NATIONAL BUREAU OF STANDARDS
microcopy resolution test chart
DIRECT AND ITERATIVE METHODS
FOR BLOCK TRIDIAGONAL LINEAR SYSTEMS
Don Eric Heller
April 1977
DEPARTMENT
of
COMPUTER SCIENCE
D D C
?f7ar?mp|z
NOV 22 19TT
gSEDTTE
Carnegie -Mel Ion University
[ON STATEMENT A
proved for public release;
Distribution Unlimited
ABSTRACT
Block tridiagonal systems of linear equations occur frequently in
scientific computations, often forming the core of more complicated prob-
lems. Numerical methods for solution of such systems are studied with
emphasis on efficient methods for a vector computer. A convergence theory
for direct methods under conditions of block diagonal dominance is developed,
demonstrating stability, convergence and approximation properties of direct
methods. Block elimination (LU factorization) is linear, cyclic odd-even
reduction is quadratic, and higher-order methods exist. The odd-even
methods are variations of the quadratic Newton iteration for the inverse
matrix, and are the only quadratic methods within a certain reasonable
class of algorithms. Serai-direct methods based on the quadratic conver-
gence of odd-even reduction prove useful in combination with linear itera-
tions for an approximate solution. An execution time analysis for a pipe-
line computer is given, with attention to storage requirements and the
effect of machine constraints on vector operations.
PREFACE
It is with many thanks that I acknowledge Professor J. F. Traub, for
my introduction to numerical mathematics and parallel computation, and for
his support and encouragement. Enlightening conversations with G. J. Fix,
S. H. Fuller, H. T. Rung and D. K. Stevenson have contributed to some of
the ideas presented in the thesis. A number of early results were obtained
while in residence at ICASE , NASA Langley Research Center; J. M. Ortega and
R. G. Voigt were most helpful in clarifying these results. The thesis was
completed while at the Pennsylvania State University, and P. C. Fischer
must be thanked for his forbearance.
Of course, the most important contribution has been that of my wife,
Molly, who has put up with more than her fair share.
ii
CONTENTS
page
1. Introduction . 1
A. Summary of Main Results 5
B. Notation 8
2. Some Initial Remarks 10
A. Analytic Tools 10
B. Models of Parallel Computation 14
3 . Linear Methods 18
A. The LU Factorization 18
1. Block Elimination 18
2. Further Properties of Block Elimination 25
3. Back Substitution 30
4. Special Techniques 31
B. Gauss-Jordan Elimination 33
C. Parallel Computation 36
4. Quadratic Methods 37
A. Odd- Even Elimination 37
B. Variants of Odd-Even Elimination 45
C. Odd- Even Reduction 48
1. The Basic Algorithm 48
2. Relation to Block Elimination 51
3. Back Substitution 56
4. Storage Requirements 58
5. Special Techniques 59
D. Parallel Computation 60
E. Summary 64
5. Unification: Iteration and Fill-In 65
iii
\
6. Higher- Order Methods 71
A. p-Fold Reduction 71
B. Convergence Rates 75
C. Back Substitution 80
7. Semidirect Methods 83
A. Incomplete Elimination 84
B. Incomplete Reduction 86
C. Applicability of the Methods 90
8. Iterative Methods 92
A. Use of Elimination 92
B. Splittings 93
C. Multi- line Orderings 94
D. Spiral Orderings 97
E. Parallel Gauss 101
9. Applications 105
A. Curve Fitting 105
B. Finite Elements 106
10. Implementation of Algorithms 110
A. Choice of Algorithms (Generalities). 110
B. Storage Requirements 113
C. Comparative Timing Analysis 120
Appendix A. Vector Computer Instructions 131
Appendix B. Summary of Algorithms and Timings 136
References 156
iv
1. INTRODUCTION
Block tridiagonal matrices are a special class of matrices which arise
in a variety of scientific and engineering computations, typically in the
numerical solution of differential equations. For now it is sufficient to
say that the matrix looks like
where b. is an n. X n. matrix and a., c. are dimensioned conformally. Thus
1 11 -N1 1 N
the full matrix A has dimension (, n.) x ( n.).
j-1 J ‘">1 J
In this thesis we study numerical methods for solution of a block tri-
diagonal linear system Ax = v, with emphasis on efficient methods for a
vector-oriented parallel computer (e.g., CDC STAR- 100, Illiac IV). Our
analysis primarily concerns numerical properties of the algorithms, with
discussion of their inherent parallelism and areas of application. For
this reason many of our results also apply to standard sequential algo-
rithms. A unifying analytic theme is that a wide variety of direct and
iterative methods can be viewed as special cases of a general matrix itera-
tion, and we make considerable use of a convergence theory for direct
methods. Algorithms are compared using execution time estimates for a
simplified model of the CDC STAR- 100 and recommendations are made on this
basis .
2
There are two important subclasses of block tridiagonal matrices,
depending on whether the blocks are small and dense or large and sparse.
A computational method to solve the linear system Ax = v should take
into account the internal structure of the blocks in order to obtain
storage economy and a low operation count. Moreover, there are impor-
tant applications in which an approximate solution is satisfactory.
Such considerations often motivate iterative methods for large sparse
systems, especially when a direct method would be far more expensive.
Since 1965, however, there has been an increased interest in the
direct solution of special systems derived from two-dimensional elliptic
partial differential equations. With a standard five-point finite differ-
ence approximation on a rectangular grid, the a^, c^ blocks are diagonal
and the b^ blocks are tridiagonal, so A possesses a very regular sparse
structure. By further specializing the differential equation more struc-
ture will be available. Indeed, the seminal paper for many of the impor-
tant developments dealt with the simplest possible elliptic equation. This
was Hockney's work [H5], based on an original suggestion by Golub, on the
use of Fourier transforms and the cyclic reduction algorithm for the solu-
tion of Poisson's equation on a rectangle. For an n X n square grid Hockney
3
was able to solve Ax * v in 0(n ) arithmetic operations, while ordinary band-
4 3
matrix methods required 0(n ) operations. Other 0(n ) methods for Poisson's
equation then extant had considerably larger asymptotic constants, and the
new method proved to be significantly less time-consuming in practice. Sub-
sequent discovery of the Fast Fourier Transform and Buneman's stable version
2
of cyclic reduction created fast (0(n log n) operations) and accurate
methods that attracted much attention from applications programmers and
r
numerical analysts [B20], [D2], [H6]. The Buneman algorithm has since
been extended to Poisson1 s equation on certain nonrectangular regions
and to general separable elliptic equations on a rectangle [S6], and well-
tested Fortran subroutine packages are available (e.g. , [S8]). Other
recent work on these problems includes the Fourier-Toeplitz methods fF2 ] ,
and Bank's generalized marching algorithms [B5]. The latter methods work
by controlling a fast unstable method, and great care must be taken to
maintain stability.
Despite the success of fast direct methods for specialized problems
there is still no completely satisfactory direct method that applies to a
general nonseparable elliptic equation. The best direct method available
for this problem is George's nested dissection [Gl], which is theoreti-
cally attractive but apparently hard to implement. Some interesting tech-
niques to improve this situation are discussed by Eisenstat, Schultz and
Sherman [El] and Rose and Whitten [R5] among others. Nevertheless, in
present practice the nonseparable case is still frequently solved by an
iteration, and excellent methods based on direct solution of the separable
case have appeared (e.g., [C4], [C5]). The multi-level adaptive techniques
recently analyzed and considerably developed by Brandt [B13] are a rather
different approach, successfully combining iteration and a sequence of dis-
cretization grids. Brandt discusses several methods suitable for parallel
computation and although his work has not been considered here, it will
undoubtedly play an important role in future studies. The state-of-the-art
for solution of partial differential equations on vector computers is sum-
marized by Ortega and Voigt [01].
References include [B19], [Dl], [S5], [S7], [S9].
4
1
With the application area of partial differential equations in mind,
the major research goal of this thesis is a general analysis of direct
block elimination methods, regardless of the sparsity of the blocks. As
a practical measure, however, the direct methods are best restricted to
cases where the blocks are small enough to be stored explicitly as full
matrices, or special enough to be handled by some storage-efficient im-
plicit representation. Proposed iterative methods for nonseparable ellip-
tic equations will be based on the ability to solve such systems efficiently.
Our investigation of direct methods includes the block LU factoriza-
tion, the cyclic reduction algorithm (called odd-even reduction here) and
Sweet's generalized cyclic reduction algorithms [S10], along with varia-
tions on these basic themes. The methods are most easily described as a
sequence of matrix transformations applied to A in order to reduce it to
block diagonal form. We show that, under conditions of block diagonal
dominance, matrix norms describing the off-diagonal blocks relative to
the diagonal blocks will not increase, decrease quadratically , or decrease
at higher rates, for each basic algorithm, respectively. Thus the algo-
rithms are naturally classified by their convergence properties. For the
cyclic reduction algorithms, quadratic convergence leads to the defini-
tion of useful semi-direct methods to approximate the solution of Ax = v.
Some of our results in the area have previously appeared in [H2].
All of the algorithms discussed here can be run, with varying degrees
of efficiency, on a parallel computer that supports vector operations. By
this we mean an array processor, such as Illiac IV, or a pipeline processor,
such as the CDC STAR- 100 or the Texas Instruments ASC. These machines
fall within the single instruction stream-multiple data stream classifica-
tion of parallel computers, and have sufficiently many common features to
L
5
make a general comparison of algorithms pc:sible. However, there are
also many special restrictions or capabilities that must be considered
when designing programs for a particular machine. Our comparisons for
a simplified model of STAR indicate how similar comparisons should be
made for other machines. For a general survey of parallel methods in
linear algebra plus background material, [H3] is recommended.
In almost every case, a solution method tailored to a particular
problem on a particular computer can be expected to be superior to a
general method naively applied to that problem and machine. However,
the costs of tailoring can be quite high. Observations about the behav-
ior of the general method can be used to guide the tailoring for a later
solution, and we believe that the results presented here will allow us
to move in that direction.
6
B[A] in the analysis of direct and semidirect methods, though it has long
been used to analyze iterative methods.
Under the assumption j| B[A]||< 1, when a direct method generates a
sequence of matrices with || B[Mi+1]|| s || B[M. ]||, we classify the method
as linear, and if || B[M^+^]|| £ || B[M^]“|| or j| B [M_^ ] 1 1 we say that it is
quadratic. The ‘“-norm is used throughout; extension to other norms would
be desirable but proves quite difficult.
The block LU factorization and block Gauss-Jordan elimination are
shown to be linear methods, and their properties are investigated in Sec-
tion 3. Sample results are stability of the block factorization when
|| B [ A ] j j < 1, bounds on the condition numbers of the pivotal blocks, and
convergence to zero of portions of the matrix in the Gauss-Jordan algorithm.
The methods of odd-even elimination and odd-even reduction, which are
inherently parallel algorithms closely related to the fast Poisson solvers,
are shown to be quadratic, and are discussed in Section 4. Odd-even elim-
ination is seen to be a variation of the quadratic Newton iteration for
A \ and odd-even reduction can be considered as a special case of Hageman
and Varga's general cyclic reduction technique [HI], as equivalent to
block elimination on a permuted system, and as a compact form of odd-even
elimination. We also discuss certain aspects of the fast Poisson solvers
as variations of the odd-even methods.
In Section 5 we show that the odd-even methods are the only quadratic
methods within a certain reasonable class of algorithms. This class in-
cludes all the previous methods as particular cases of a general matrix
iteration. The relationship between block fill-in and quadratic conver-
gence is also discussed.
L
7
In Section 6 the quadratic methods are shown to be special cases of
families of higher-order methods, characterized by || B[M^+^]|| ^ j| BTM^ ] 1 1 P-
We note that the p-fold reduction (Sweet's generalized reduction [S10]) is
equivalent to a 2- fold reduction (odd-even) using a coarser partitioning
of A.
Semidirect methods based on the quadratic properties of the odd-even
algorithms are discussed in Section 7. Briefly, a semidirect method is a
direct method that is stopped before completion, producing an approximate
solution. We analyze the approximation error and its propagation vis-a-vis
rounding error in the back substitution of odd-even reduction. Practical
limitations of the semidirect methods are also considered.
In Section 8 we consider some iterative methods for systems arising
from differential equations, emphasizing their utility for parallel com-
puters. Fast Poisson solvers and general algorithms for systems with
small blocks form the basis for iterations using the natural, multi-line
and spiral orderings of a rectangular grid. We also remark on the use of
elimination and the applicability of the Parallel Gauss iteration [T2],
[H4 ] .
A brief section on applications deals with two situations in which
the assumption |(B[A]j|< 1 is not met, namely curve fitting and finite ele-
ments, and in the latter case we describe a modification to the system
Ax = v that will obtain ||b[A]|| < 1.
Methods for solution of block tridiagonal systems are compared in Sec-
tion 10 using a simplified model of a pipeline (vector) computer. We con-
clude that a combination of odd-even reduction and LU factorization is a
powerful direct method, and that a combination of semidirect and iterative
methods may provide limited savings if an approximate solution is desired.
8
A general discussion of storage requirements and the effect of machine
constraints on vector operations is included; these are essential con-
siderations for practical use of the algorithms.
l.B. Notation
Let the block tridiagonal matrix
Here b. is an n. X n. matrix and a. = 0, c = 0. The subscript N in the
i l l I N
triplex notation will often be omitted for simplicity. Define
L[A] = (a., 0, 0) ,
D[ A ] = (0, b , 0),
U[A] - (0, 0, c ),
B[A] = I - (D[A])_1A = (-b"1 a , 0, -b"1 c ) ,
C[A] - I - A(D[A])” 1 = (-a b“}i, 0, -c jbj^) .
Use of these notations for general partitioned matrices should be evident;
nonsingularity of D[A] is assumed in the definitions of B and C. B[A] is
the matrix associated with the block Jacobi linear iteration for the solu-
tion of Ax = v. It also assumes an important role in the study of certain
semidirect methods and will be central to our analysis of direct methods.
9
We further define = diag ( 0 , . . . , 0 , I, 0,...,0), where the I is
Oj X n^ and occurs in the jth block diagonal position. Pre- (post-) mul-
tiplication by E. selects the jth block row (column) of a matrix, so that,
J N N
E . AE . , D[A]E. = E.U[A], and E. L[A ] E . = 0.
. , l i j 1 t,. . i j
i= 1 J J j=i J
for example, D[A] =
For an nXn matrix M and an n-vector v, let
Pi(M)
y. (M)
II M|| =
II Ml'
n
"j-i
n
'i*l
m. .
1 1
|m. . = p . (M ) ,
J-J1 J
maxi * c^e co"nonn>
max ^ v (M)
M
v|| = max. |v.J ,
n
i=l
Vi'
|M
( I mi j I ) » and
0 (M) = the spectral radius of M.
The reader should take care not to confuse row sums (p^(M)) with the spec-
tral radius (o(M)). When A = (a., b., c .) , let \(A) = max . (4 ||b . ^ a . || c-
111 j J J J J
The condition \(A) £ 1 will be used as an alternative to |j B[A]|| < 1.
-1
2. SOME INITIAL REMARKS
2 . A . Analytic Tools
Much of the analysis presented in later sections involves the matrix
B[A]. Why should we consider this matrix at all?
First, it is a reasonably invariant quantity which measures the diag-
onal dominance of A. That is, if E is any nonsingular block diagonal matrix
conforming with the partitioning of A, let A* 3 EA , so B* = B[A*] = I-D* *A*
= I - (ED) (EA) = I - D ^A = B. Thus B[A] is invariant under block scal-
ings of the equations Ax = v. On the other hand, B is not invariant under
a change of variables y = Ex, which induces A' = AE We have
B' = B[A' ] = EBE \ and for general A and E only the spectral radius 3(B)
is invariant. Similarly, C[A] measures the diagonal dominance of A by
columns, and is invariant under block diagonal change of variables but not
under block scaling of the equations.
Secondly, one simple class of parallel semidirect methods (Section 7)
can be created by transforming Ax = v into A*x = v*, A* = MA, v* = Mv for
some matrix M, where || B[A*]|| « || B[A]||. This transformation typically
represents the initial stages of some direct method. An approximate solu-
tion y = D[A*] *v* is then computed with error x - y = B[A*]x, so ||B[A*]||
measures the relative error. Thus we study the behavior of B[A] under
various matrix transformations.
Thirdly, B[A] is the matrix which determines the rate of convergence
of the block Jacobi iteration for Ax = v. This is a natural parallel itera-
tion, as are the faster semi- iterative methods based on it when A is posi-
tive definite. It is known that various matrix transformations can affect
the convergence rates of an iteration, so we study this property.
11
Finally, B[A] and C[A] arise naturally in the analysis of block elim-
ination methods to solve Ax = v. C represents multipliers used in elimina-
tion, and B represents similar quantities used in back substitution. In
order to place upper bounds on the size of these quantities we consider
the effects of elimination on || B[A]||.
Since we will make frequent use of the condition jj B[A]|| < 1 we now
discuss some conditions on A which will guarantee that it holds.
Lemma 2.1. Let J = + J? be order n, |j| = | J^j + | J? | , and suppose K
satisfies K = J^K + J9 . Then || j|| < 1 implies p^CK) £ P^J) for 1 £ i £ n,
and so || K || £ j| j||.
Proof . From |j^| + j J0 1 = |j| we have p^J^) + P^C^) = Pj_(J) for 1 £ i £ n.
From K = J^K + J we have p^(K) £ o (J^) || K | j + oi(J?). Suppose |j k|| s 1.
Then P^K) £ C ^ ( JT^) || k|| + P.^^) II Kll = II K|| s || J|| || K||> which implies
|| K || £ || J j| j| K jj and 1 £ || j||, a contradiction. Thus j| k|| < 1 and
Pl(K) £ O.^) + p^J.,) = Pt(J). QED.
We note that if o ^ ( J * 0 then we actually have p^K) < p^(J) and
if p.(J,) = 0 then p . (K) = o.(J„) = p,(J). It follows that if each row of
l 1 iv t 2 i
contains a nonzero element then || K|| < || j||. These conditions are rem-
iniscent of the Ste in-Rosenberg Theorem ([V3], [H8]) which states that if
J is irreducible, and nonnegative, J9 * 0 and p(J^) < 1 then p(J) < 1
implies p (K) < p(J) and p(J) = 1 implies p(K) = o(J). There are also cor-
respondences to the theory of regular splittings [V3].
Lemma 2.2. Suppose J = (I - S)K + T is order n, |j| = | (I - S) K | + |t|,
P-^S) £ P^(T) , 1 £ i £ n. Then || j|| < 1 implies p^K) £ p^J) , 1 £ i £ n,
and so f| K|| £ || j|j.
12
Proof. For 1 £ i £ n, p^T) £ p.(T) + p.((I - S)K) = P^J) £ || j|| < 1,
so p.(T) < 1. Now, 1 > 0i(J) = p.((I - S) K) + Pt(T) S pt(K) - p.(S)||K||
+ p.(T) & P.(K) + P.(T)(1 - || K|| ). Suppose || K|| = p^(K) > 1. Then
I > pj(K) + p (T)(l - pj(K)), and we have 1 < p^(T), a contradiction.
Thus || K || ^ 1 and pi(J) 2 Pi(K). QED.
If S - T in Lemma 2.2 then K=»SK+ (J- S), |j| = | (I - S)K|+ |s|
= | J - S j + | S | , so by Lemma 2 . 1 we again conclude || K|| £ || j||.
Similar results hold for the 1-norm.
Lemma 2.3. Let J * J + J,, , |j| 3 |j | + | J0 | , and suppose K satisfies
K - KJX + J . Then || J 1^ < 1 implies || K||x ^ || j||x .
Proof. Since JT = + J^, |jT| => | J* | + | J* | , KT = J^KT + J, and
II Ml^ = || MT||, we have || JT|| = || < 1, and by Lemma 2.1, || - || KT||
* (I JT|( " (1 J|lr QED-
Lemma 2.4. Suppose J = K(I - S) + T, | j| = | K(T - S) | + |t| , v^S) £ y^T) ,
1 * i £ n. Then || j||. < 1 implies || k|L * || j|| .
Proof. We have v^M) = P^fM ), = (I - ST)KT + TT, so by Lemma 2.2
|| K^|| £ || J^||, and the result follows. QED.
To investigate the condition || B[A]||< 1, we first note that if A is
strictly diagonally dominant by rows then j| B[A]|| < 1 under any partition-
ing of A. Indeed, if we take J to be the point Jacobi linear iteration
matrix for A, J 3 I - diag(A) ^A, * D[J], Jj * J “ then
|| J1|| * || J|| < 1, |J| = | JL | + | J2| and B[A] - (I - so
|| B[A]|| * i| J|| < l by use of Lemma 2.1.
Now, if we denote the partitioning of A by the vector II = (n^ ,n,, , . . . ,n^)
then the above remark may be generalized. Suppose H' = (m^ ,ra9 , . . . .m^) rep-
resents another partitioning of A which refines II. That is, there are inte-
Pi_1
gers 1 = pQ < p^ < . . , < pN = M4-1 such that ni = m.. If B[A ] and
1 ~ J = P i- 1 J
B'[A] are defined by the two partitions and j| B'[A]|| < 1 then [| B [ A ] 1 1
£ || B'[A]||. This is an instance of the more general rule for iterative
methods that "more implicitness means faster convergence" (cf. [H8, p. 107])
and will be useful in Section 6. Briefly, it shows that one simple way to
decrease || B[A]|| is to use a coarser partition.
Let us consider a simple example to illustrate these ideas. Using
2
finite differences, the differential equation -u - u + Xu = f on T 0 , 1 ] ,
n xx yy l > j >
X > 0, with Dirichlet boundary conditions, gives rise to the matrix
9
A = (-1, M, -I), M = (-1, 4+Xh", -1). By Lemma 2.1, || B[A]|| £ || J[A]|| =
2
4/ (4+Xh ) < 1. Actually, since B[A] is easily determined we can give the
better estimate || B[A]|| = 2 1 1 M ^ || < 2/ (2+Xh~) .
Indeed, many finite difference approximations to elliptic boundary
value problems will produce matrices A such that || B[A]|| < 1 [ V3 ] ; |{ B[A]||
will frequently be quite close to 1. This condition need not hold for
certain finite element approximations, so a different analysis will be
needed in those cases. Varah [Vl], [V2] has examples of this type of
analysis which are closely related to the original continuous problem. In
Section 9.B we consider another approach based on a change of variable.
|| B[ A ] 1 1 < l is a weaker condition than strict block diagonal dominance
relative to the “-norm, which for block tridiagonal matrices is the condi-
tion ([FI], [Vl]) || b j 1 1| (|| a j || + || c j 1 1 ) < 1, 1 * j s u. The general
definition of block diagonal dominance allows the use of different norms.
so for completeness we should also consider theorems based on this
14
assumption. However, in order to keep the discussion reasonably concise
this will not always be done. It is sometimes stated that block diag-
onal dominance is the right assumption to use when analyzing block elimi-
nation, but we have found it more convenient and often more general to
deal with the assumption || B[A]||< 1.
Other concepts related to diagonal dominance are G- and H-matrices .
If J = I - diag(A) ^A= L + U, where L (U) is strictly lower (upper) trian-
gular, then A is strictly diagonally dominant if || j|| < 1, a G-matrix if
|| (I - j L | ) 1 |U| || < 1, and an H-matrix if p(|j|) < 1. Ostrowski [02]
shows that A is an H-matrix iff there exists a nonsingular diagonal matrix
E such that AE is strictly diagonally dominant. Varga [V4] summarizes
recent results on H-matrices. Robert [R3] generalizes these definitions
to block diagonally dominant, block G- and block H-matrices using vector-
ial and matricial norms, and proves convergence of the classical iterative
methods for Ax = v. Again, for simplicity we will not be completely gen-
eral and will only consider some representative results.
2 .B . Models of Parallel Computation
In order to read Sections 3-9 only a very loose conception of a vector-
oriented parallel computer is necessary. A more detailed description is
needed for the execution time comparison of methods in Section 10. This
subsection contains sufficient information for the initial part of the
text, while additional details may be found in Appendix A, [H3], [S3], and
the references given below.
Parallel computers fall into several general classes, of which we con-
sider the vector computers, namely those operating in a synchronous manner
and capable of performing efficient floating point vector operations.
15
Examples are the array processor Illiac IV, with 64 parallel processing
elements [B6], [B12]; the Cray-1, with eight 64-word vector registers [C8];
and the pipeline processors CDC STAR- 100 [C2] and Texas Instruments ASC
[Wl], which can handle vectors of (essentially) arbitrary length. These
machines have become important tools for large scientific computations,
and there is much interest in algorithms able to make good use of special
machine capabilities. Details of operation vary widely, but there are suf-
ficiently many common features to make a general comparison of algorithms
poss ible .
An important parallel computer not being considered here is the asyn-
chronous Carnegie-Mellon University C.mmp [W4], with up to 16 minicomputer
processors. Baudet [B7] discusses the theory of asynchronous iterations
for linear and nonlinear systems of equations, plus results of tests on
C . mmp .
Instructions for our model of a vector computer consist of the usual
scalar operations and a special set of vector operations, some of which
will be described shortly. We shall consider a conceptual vector to be
simply an ordered set of storage locations, with no particular regard to
the details of storage. A machine vector is a more specialized set of
locations valid for use in a vector instruction, usually consisting of
locations in arithmetic progression. The CDC STAR restricts the progres-
sion to have increment = 1, which often forces programs to perform extra
data manipulation in order to organize conceptual vectors into machine
vectors .
We denote vector operations by a parenthesized list of indices. Thus
the vector sum w = x+y is
16
w. = x. + y . , (1 s i s N)
1 t i
, N
if the vectors are in R , the componentwise product would be
Wi = Xi X yi’ C1 s ^ N> »
and a matrix multiplication C = AB would be
Cij=_^1aik v (1 "N; 1 sn)-
Execution time for a vector operation depends on the length of the vector
and the internal mechanism of the vector computer, but can usually be
written in the form T N + c if the vector length is N. On a pipeline
op op r r
computer, 1/ t is the asymptotic result rate and a is the vector instruc-
op op
tion startup cost. On a k-parallel computer (one with either k parallel
processors or k-word vector registers) the basic time would be t for k
op
operations in parallel and the vector instruction cost t N k~ + overhead.
op
In either case the tN + a model is adequate for vector operations; scalar
operation times are given as t
op
Simplified instruction execution times
for the CDC STAR- 100 are given below, in terms of the 40 nanosecond cycle
time. Usually T « t « a
op op op
operation
scalar time
vector
time
add, subtract
13
b*
71
multiply
17
N +
159
divide
46
2N +
167
square root
74
2N +
155
transmit
11
91
maximum
-
6N +
95
summation
-
6.5N
+ 122
inner product
-
6N +
130
17
The use 01 .ong vectors is essential for effective use of a vector
computer, in order to take full advantage of increased result rates in
the vector operations. The rather large vector startup costs create a
severe penalty for operations on short vectors, and one must be careful
to design algorithms which keep startup costs to a minimum. Scalar opera-
tions can also be a source of inefficiency even when most of the code con-
sists of vector operations. Data manipulation is often an essential con-
sideration, especially when there are strong restrictions on machine vec-
tors. These issues and others are illustrated in our comparison of algo-
rithms in Section 10, which is based on information for the CDC STAR- 100.
18
3. LINEAR METHODS
The simplest approach to direct solution of a block tridiagonal linear
system Ax = v is block elimination, which effectively computes an LU factor-
ization of A. We discuss the process in terms of block elimination on a
general partitioned matrix under the assumption j| B[A]|j< 1; the process is
shown to be linear and stable. Various properties of block elimination are
considered, such as pivoting requirements, bounds on the condition numbers
of the pivotal blocks, error propagation in the back substitution, and
special techniques for special systems. Gauss- Jordan elimination is studied
for its numerical properties, which foreshadow our analysis of the more impor-
tant odd-even elimination and reduction algorithms.
3 .A. The LU Factorization
3.A.I. Block Elimination
The block tridiagonal LU factorization
A = (aj,b.,c.)N = (i.,I,0)N(0,d.,c.)N
is computed by
d = b
1 1
,-l
dj = bj - ajdj-l cj-l> J = 2.....N,
l3 = Vj-l’ j =
where the existence of d^. ^ is assumed and will be demonstrated when || B [ A ] 1 1 < 1.
For the moment we are not concerned with the particular methods used to compute
or avoid computing either Z ^ or d^ ^ cj-l"
19
Let A^ = A and consider the matrix A0 which is obtained after the first
block elimination step. We have A0 = (a.,3.,cJN> 3^ = d^ = b^, a ^ = 0,
32 = d2’ aj = aj and = bj for 3 s Let Bx = B[A1], B2 = B[A,].
Note that and differ only in the second block row; more specifically,
in the (2,1) and (2,3) positions.
Theorem 3.1. If J| B ^ || < 1 then || B2 || £ || B || .
Proof. Let T = B^E^ and let S = TB^. We have d2 = b9(I - b0^a0b^c ),
II b2 la2bllcl II = H SH ^ II B1 II II Ei II II Bill = II Bil|2 < 1> so d21 exists and B„
is well-defined. Also, B^ = (I - S)B, + T, { B^ | = j (I - S) B, | + j T | by con-
struction, and o^(S) £ 0 j (T) 1 1 B^|[ ^ P^(T). Thus Lemma 2.2 applies, yielding
QED
Now, for any partitioned matrix the N stages of block elimination (no
block pivoting) can be described by the process
for i = 1 N- 1
Ai+1 - (I ' L[Ai]EiD[A1]‘1)Ai
= (I + L[C[Ai]]Ei)A..
We transform A^ into by eliminating all the blocks of column i below the
diagonal; when A is block tridiagonal A^ is the upper block triangular matrix
(0,dj,Cj)^. The same effect can be had by eliminating the blocks individually
by the process
20
Aj = A
for i = 1, . . . ,N- 1
Af1! = A
l+l l
for j = i+l,...,N
A^j) = (I - E.Af-^E.DCA^"^ ]'1)A^"1)
l+l j i+l l l+l J l+l
L = (I + EjC[Ai+r ]Ei
i+l
L
A . . = a<N>
i+l i+l
Each minor stage affects only block row j of the matrix. It is seen that
block row j of A^~ ^ = block row j of A.,
i+l J l
block row i of Ap,1"* = block row i of A..
i+l l
from which it follows that
E.A^YofAy- Y1 = E.A.E.D[A.]-1,
j i+l l i+l J j l l L iJ ’
so the major stages indeed generate the same sequence of matrices A^. Let
B. = B[A.].
Corollary 3.1, If A is block tridiagonal and j| B^|| < 1 then || B.+1|| ^ |j B. ||,
1 £ i SN-1,
The proof of Corollary 3.1 requires a bit more care than the proof of Theorem
3.1 would seem to indicate; this difficulty will be resolved in Theorem 3.2.
We should also note the trivial "error relations"
*i+i - 'Hill s ll*i - *JI.
Bi+1- 'nII sliBi- BJ-
21
These follow from the fact that, for block tridiagonal matrices,
i N
A.
1
E .
A
+
~j=i+l
EjAr
, which holds regardless of the value of j| B ^|| .
Wilkinson [W3] has given a result for Gaussian elimination that is similar
to Corollary 3.1, namely that if A is strictly diagonally dominant by columns
then the same is true for all reduced matrices and || A^+^ (reduced) ||^ £
|| A^ (reduced) j|^. This actually applies to arbitrary matrices and not just
the block tridiagonal case. Brenner [B14] anticipated the result but in
a different context.
Our next two theorems, generalizing Corollary 3.1 and Wilkinson's observa-
tion, lead to the characterization of block elimination as a norm- non- increas ing
k
projection method , and proves by bounding the growth factors that block
elimination is stable when ))b[A]|| < 1. In addition, many special schemes for
block tridiagonal matrices are equivalent to block elimination on a permuted
matrix, so the result will be useful in some immediate applications.
Theorem 3.2. If )| B 1|| < 1 then || B.+1|| <; || B.|] and j| At+1|| ^ || A.|| for
1 £ i S8-1.
Proof. Suppose || B^||< 1. We have
A1+1 3 (I - L[A.]EiD[Aif1)Ai
- A. - L[A.]E. + L[A1]E1B1,
PJ(Ai+l>* Pj(Ai ' + Pj(L[Ai]EiB.)
s P j(Ai - L[A.]E.) + oj(L[Ai]Ei)|| B.||
k
The projections are onto spaces of matrices with zeros in appropriate positions.
* Pj(A. - L[A.]E.) + PjCLCA.JE.)
= P.(A.),
S° II Ai+lH * II Ai II • Now let
Q = D[Ai]‘1L[Ai]EiD[Ai]'1A.>
30 Ai+1 = Ai ' °fAi^» D[Ai+l] = D^AiKI - D[Q ]) •
Let S = D[ Q ] - Since
Q = (-L[B.]E.)(I - B.) = -L[B.] + L[B.]E.B.,
S = D[L[B.]E.B.] and ||s||* ||L[B.]|| ||e.|| ||B.||s ||B.j|2 < 1. Thus
D[A^+^] exists and B^+^ is well-defined. It may be verified that
(I ' = Bi - S + Q. Let T = Bl - S + Q, J = S + T = B + Q. Since
D[ Bjl ] = D[Q - S] = 0, we have D[T] = 0 and |j| = |s| + |t|. Thus, by Lemma 2.1,
j| j|| < 1 implies || Bi+1|| ^ j| j|| . But we have, after some algebraic manipula-
tion.
J = B. + Q
i
Ci EA
yN
+ / E.B, (I - E. + E.B.) ,
Tc-i+1 k 1 1 11
and it follows from this expression that || j|| £ || B^|| . QED
Although expressed in terms of the N- 1 major stages of block elimination,
it is clear that we also have, for B^ = B[A^ ], || B^ || £ || B^"1^ || and
llAij)|| 5 || X) ||. Moreover, as long as blocks below the diagonal are
eliminated such inequalities will hold regardless of the order of elimination.
23
We shall adopt the name "linear method" for any process which generates
a sequence of matrices such that || B[M^+^][| s || B[M^]||, given that
|| B[M^]||< 1. Thus block elimination is a linear method.
A result that is dual to Theorem 3.2 is the following.
Definition. Let CL~ be the reduced matrix consisting of the last N-i+1 block
rows and columns of A^. = A^ = A) .
Theorem 3.3. If || CfA]^ < 1 then || C[<2.+1 ] || l * || C[£2. ] ||r \\^i+l\\l ^ H^llj.
and || L[C1]E.||1 * || CCA]^, C. = C[A.].
Remark. It is not true in general that ||c || £ | j c ^ 1 1 f HCilli< 1’
Proof. Let the blocks of C . be c , 1 £ ra, p £ N. Then C\&. 1 = (c ),-<s„
i mp L mp
mp
Let K = (cmp + craicip) i+l^n,p^
= ( r ) , , •+■
' mp i+ls<n,p^I
(ci,i+l CiN)-
Using Lemma 2 .4 ,
But
^i+iHl!
if : k L < l.
VK>*V E C6J? ])
4 H ~k=i+l
i+l»i\
'Ni
l^qCEf C82,])
+ II cBS1]||1vq(EiCK?i])
s v-C1+i E*
r
24
1
+ vq(Et CE^])
= vq(C52.]).
Thus II k||x ^ ||c[<2i]||1< 1 and || ct^i+1] II ± 5 || C^] || 1 . By construction
this implies || £ || C [ A ] 1 1 ^ for 1 £ i £ N. Now, let the blocks of
be a , 1 ^ m, p s n. Then CL. , = (a + c .a. and
mp l+l mp mi ip i+l^n,p^l
s V^+1 Etf>
llVq(Ei°l)
' X, «V ■
QED
We note that if A is symmetric then C[A] = B[A ] and Theorems 3.2 and 3.3
become nearly identical. Of course, B and C are always similar since C =* DBD
When A is positive definite, Cline [ C 3 ] has shown that the eigenvalues of 42
interlace those of 42^ Actually, the proof depends on use of the Cholesky
factorization, but the reduced matrices are the same in both methods. In
consequence of Cline's result we have ||^Z || 5 •
The statements of Theorems 3.2 and 3.3 depend heavily on the particular
norms that were used. It is not true in general, for example, that if
||B[A1]||1< 1 then || B[A2] ||j £ || B[Aj] ||j or even || B[£, ] ^ s || B[4^ ] ||r This
makes the extension to other norms quite difficult to determine.
25
1
3. A. 2. Further Properties of Block Elimination
It is well-known that the block factorization A - LAU is given by
N
L - I - ; L[c ]E ,
S-i J J
A = D[A^1,
N
U = I - B =* I - / E. U[B ].
N “j-l J J
We have therefore shown in Theorems 3.2 and 3.3 that if || 1 then
II L||x * 1 + II C1|| J_ <2, || L|| * 1 + Nn|| C1||1, n = maXj n . , and if || B1|| < 1
then || U|| * 1 + || B x j| <2, || u|| l £ 1 + Nn|| bJI . (For block tridiagonal
matrices the factor N can be removed, so |j L|| £ 1 + nj| C ||^, || £ 1 + n|| B^j|.)
Moreover, the factorization can be completed without permuting the block rows
of A. The purpose of block partial pivoting for arbitrary matrices is to
force the inequality || L[ C ^ ]E ^ 1 1 £ 1. Broyden [B15] contains a number of
interesting results related to pivoting in Gaussian elimination. In particular
it is shown that partial pivoting has the effect of bounding the condition
number of L, though the bound and condition number can actually be quite
large in the general case. We will have some more comments on pivoting after
looking at some matrix properties that are invariant under block elimination
(cf. [311]).
Corollary 3.2. Block elimination preserves strict diagonal dominance.
Proof . It follows from Theorem 3.2 that Gaussian elimination preserves strict
diagonal dominance: simply partition A into lxl blocks. Now suppose we have
A^ and are about to generate A This can be done either by one stage of
block elimination or by n^ stages of Gaussian elimination, for both processes
produce the same result when using exact arithmetic [HS, p. 130]. Thus if A^
is strictly diagonally dominant, then so is QED
r
26
It is shown in [B4] that Gaussian elimination maps G-matrices into
G-matrices. This can be extended in an obvious fashion to the case of block
elimination. We also have
Corollary 3.3. Block elimination maps H-matrices into H-matrices.
Proof. Recall that A is an H-matrix if and only if there exists a nonsingular
diagonal matrix E such that A' = AE is strictly diagonally dominant. Con-
sider the sequence
A^ = A' ,
A^+1 = (I - L[A'.]E. D[A’]-1)A^.
Clearly A^ = AjE; suppose that A^ = A^E is strictly diagonally dominant. We
then have
Ai+1 = (I " L[Ai]E EiE“1D[Aif1)A1E
= (I - L[At]E. DrA.f^A.E
= Ai+1E’
so A^+^ is an H-matrix since A^_^ is strictly diagonally dominant by Corollary
3.2. QED
We have chosen to give this more concrete proof of Corollary 3.3 because
it allows us to estimate another norm of B[A] when A is an H-matrix. Define
I|m|!e= || E ^ME 1 1. Then since A' = AE implies B[A' ] = E ^B[A]E it follows
that if A is an H-matrix then || BCAi+1]||E * || BCA^H . Since the matrix E
will generally remain unknown this is in the nature of an existential proof.
The proof of Corollary 3.2 brings out an important interpretation of
Theorem 3.2. If block elimination is actually performed then our result
gives direct information about properties of the reduced matrices. If not,
and ordinary Gaussian elimination is used instead, then Theorem 3.2 pro-
vides information about the process after each n^ steps. While we are not
assured that pivoting will not be necessary, we do know that pivot selec-
tions for these n. steps will be limited to elements within the d. block
for the block tridiagonal LU factorization. When A has been partitioned to
deal with mass storage hierarchies this will be important knowledge (see
Section 10. A).
One more consequence of Theorem 3.2 is
Corollary 3.4. In the block tridiagonal LU factorization, if Jj B^|j < 1 then
cond(dj) = || d || || d" 1 1| < 2 cond(b.) / (1 - || B1||2 ).
Proof. Letting a , = b . *a .d . \c . . , we have d. = b.(I - a.), ||a\|| £
J J J J-l J-l J J j 11 j"
II Bx|| || Bn|| * || Bjf < 1, so that || d || s || b || (1 + || o || ) < 2 || b ||, and
* b
(I- a.)"1!!* Hb-.1!!/ (1- Hcjll) ^ || b”1 1| / (1- Hb^I2).
It is important to have bounds on cond(d^) since systems involving d^
must be solved during the factorization. Of course, a posteriori bounds on
|| d.^|| may be computed at a moderate cost. We note also that
cond(A) ^ maXj 2 cond(b^) / (1 - || B^H ) when || B^|| < 1, and it is entirely
possible for the individual blocks b. and d. to be better-conditioned than A.
J j
The inequalities in Theorem 3.2 are the best possible under the condition
|| B^ll < 1. This is because, for example, the first block rows of B^ and B?
are identical. However, it is possible to derive stronger results by using
different assumptions about B.. .
V
28
i
Theorem 3 .4 . [H4 ] For A = (a^ ,b ,c ^) , let \(A) = max^ (4 || *a^ || ||b^ ^c ^ || ) ,
e j = b'V. If X(A) S 1 then || I - e j || £ (1 - Jw\)/2 and || 1 1| <: 2/(1 + v^) .
Some useful consequences of this theorem include better estimates for
II B., || and a smaller upper bound on the condition numbers of the d.'s. For we
N J
have B^ = (0,0,-dj^cJ = (0,0, -e^b^ ^c^) , which implies that
|| B^jU ^ 2 || U[B^ ] || / (1 + v/1- X) . As an example , if A = (1,3,1) then Theorem
3.2 predicts only |) B^H ^ || B^|| = 2/3, while by Theorem 3.4 we obtain
|| Bj| £ 2(1/3) / (1 + ,/l-4/9) * 0.382. Also, || d“ 1 1| £ || d~ Lb ^ || ||b~ 1 1| =
l|a'j1|| l|b j 1 II S -||b j 1 ||» and ||d. || * |^- d . || + |^ || - INjd “ -^ || + |^ || ^
ll b j || || I - e ^ |j + || b. j| £ 3 || b^l/2. It follows that cond(dj) = || d" X|1 ||d j || s
3 cond(bj) . More precisely, cond(dj) ^ (1 + 2 u ( A.) ) cond(bJ , u * ( 1 - • '1- \)/( 1 + »/l- X) .
This is further proof of the stability of the block LU factorization.
Theorem 3.4 has been used in the analysis of an iteration to compute the
block diagonal matrix (0,d^,0)^ [ H4 ] . Such methods are potentially useful
in the treatment of time- dependent problems, where good initial approximations
are available. This particular iteration is designed for use on a parallel
computer and has several remarkable properties; for details, see [H4] and
Section 8.E.
Two results on the block LU factorization for block tridiagonal matrices
were previously obtained by Varah.
Theorem 3.5. [ Vl ] If A is block diagonally dominant (|| b^||(||a^||+ | (c ^ 1 1) s 1)
and || Cj|| J* 0, then || d~ 1 1| * || c^f1, || || s || || / || and
11^11*11^11+ II aj II*
29
In comparison to Theorem 3.2, the bound on || d^|| corresponds to the
bound || < w^ich follows from || B^U < 1, and the bound on || d^ ||
corresponds to the bound ||a^|| ^ ||a||. In contrast to Theorem 3.2, j| B || < 1
“1 "1 2
only allows us to conclude || d^ || 5 || b^ || / (1 - j| B^|| ), and we can only
estimate || ^ || in the weak form || ijcj_1|| * || a. || || dj^c^JI 5 || « ||.
Varah's second result is
Theorem 3.6. [Vl] Let a. = (|| bj^ajll II ^j\Cj ill 2 £ j £ N, and sup-
pose a. P 0. Then the block LU factorization of A is numerically stable
J
(i.e., there exists a constant K(A) such that max(|| l || , || d^ || ) £ K(A)) if
(cty 1 , o' j is positive semidef inite.
We note that (<*^,1,3. is a symmetrization of ( 1 1 b . ^a ^ ||, 1 , | j b ^ ^c ^ 1 1 )
and that \(A) £ 1 implies that it is positive definite. The actual bounds
on || 2^ ||, || d j 1 1 given by Varah [VI] for this case will not be important for
this discussion.
Since many block tridiagonal systems arise from the discretization of
differential equations, a number of results have been obtained under assump-
tions closely related to the original continuous problem. Rosser [R7] shows
that if the system Ax = v, A = (aj,bj,Cj), satisfies a certain maximum
principle then A is nonsingular, the LU factorization exists (each d^ is non-
singular) , B^ is non-negative and ||b^ || S 1. Improved bounds on B^ are
obtained when A is defined by the nine-point difference approximation to
Poisson's equation. In particular there is the remarkable observation that
the spectral condition number of d^ is less than 17/3. Such results depend
strongly on the particular matrix under consideration. For examples of
related work, see [B18], [Cl], [D2].
30
3. A. 3. Back Substitution
So far we have only considered the elimination phase of the solution of
Ax = v. For block tridiagonal matrices A = (a^.b^.c^^ this consists of
the process
di = V fi = v
d. = b - a.d.^.c. f. = v. - a.d.^.f.
J j J 3-1 3-1 3 3 3 3-1 j-1*
3 = 2 , . . • ,N.
The back substitution phase consists of the process
XN = dN fN’
X3 = dj1(f3 " CjXJ+l)* j =
In matrix terms these are the processes
Al - A, f(1) = v,
A.+1 - (I + L[C.]E.)A. , i =* 1,...,N-1,
f(L+1) = (I + L[C.]E.)f(L) , i = 1,...,N-1,
(N) „r» n-lr(N>
g = x D[A^] fv ,
x(l) =» g + EjB x(l+1), i =* N- 1 1 ,
X 53 X
(1)
An alternate form of back substitution is
x(N) - D[AN]’1f(N)
.(1)
X 3 X
(I + BNE. + 1)x(i+1), i = N- 1 , . . . , 1 ,
(1)
31
In either case it is seen that plays an essential role in studying the
back-substitution process.
For those situations where a band elimination method for Ax = v might be
preferred, the matrix maintains its importance for back substitution. We
have A = (a.,b.,c.) = (a . , l . ,0) (0,u . ,c .) , d. = i.u., a. = a .u . _, c. = l, c.,
JJJ J J JJJ JJJ J J-l J JJ
A A
where .^(u^) is lower (upper) triangular. The elimination phase generates
* * 1 * - 1
vectors f. = 2.f., and the back substitution computes x^ = u^ f^ = f ,
A-l * A . -1
x. = u. (f. - c.x.,,) = d. (f. - c.x.,..). Except for the specific form of
J J V J J J+l J J J J+l
roundoff errors the band and block back substitutions must behave identically.
Let us briefly consider the effect of rounding errors in back substitu-
tion; the computed solution will be y . = x . + * and we require bounds on f .
Suppose the computed forward elimination yields d^ = d^ + 5 , f^ = f^ -4- cp ^ .
JU
Let 7). represent the errors introduced in forming f. - c.y. ,, so that y.
J J J J+l J
* (2) *
satisfies (d . + 6. )y. = f. - c.y. , + T|. if ordinary Gaussian elimination
J J J J J J+l J
*
is used to solve for y. in d. y. = f. - c.y.,, + T|.. It follows that
J J 'J J J J+l J
v. = d.^(f. - c.y.,,) + e., where e. = d.^(®. + 1. - (6^ + 6^ )y.), which
J J J J J+l j j J J j J j j
implies = - djlcj?j+l- Thus H Sjll 5 II BNH II ?j+lH+ II ejH» and the
previous errors tend to be damped out in later computations if ||b^||< 1. If
we have a uniform upper bound for the e^'s, || || £ e, and II Bj,|| < 1, then
we
obtain J! Sj|| 5 (N - j + 1) e> || J] * d - || j+1) e/ d - \^\\> < •/ (1- |)bn|J)
3. A. 4. Special Techniques
In some cases the matrix A = (a^b^c^) Possesses special properties
that may be important to preserve in order to satisfy a strong stability
criterion. For example, we might have b
j
t . - a , - c . where
J j J
t j || is small
or zero; this requires that all blocks be n X n. Scalar tridiagonal systems
J
r
32
with this property arise from the discretization of second order ODE boundary
value problems. Babuska [Bl-3] has discussed a combined LU - UL factoriza-
tion for such a matrix which allows a very effective backward error analysis
to the ODE and a stable numerical solution of a difficult singular pertur-
bation problem posed by Dorr [ D3 ] . (See [H7] for a symbolic approach to
this problem.)
There are two essential features to Babuska' s method. The first is
to express b^ and d^ indirectly by starting with t ^ and computing the se-
quence s^ = t^, Sj = tj - ]_• By this we have d = s. - c^. The
second feature is to augment the forward elimination with a backward elim-
ination and dispense with the back substitution by combining the results
of elimination. Specifically, suppose A = (— p j , t ^ + p^. + q^.-q^)^.
Pj.q^.tj £ 0, p^ = q^ = 0, A nonsingular. The algorithm for Ax = v is
V fl
V SN
V fN = V
-1
Sj = Cj + Pj(Sj-l+ qj-l} SM
£j ' vj + pj<sj-i +
★ ★ _ \ ★
s, = t + q .(s , , . + p^.) s
j
*j 'j+l rj+i “j+i
f. ■ v + q . (s . . + p. .
J j V J+l KJ+1 J+l
j 2 , . . . ,N,
j - N-
(s. + s* - tj)_1(fj + f * - v.), j = 1, . . . ,N.
In matrix terms, we factor A = A^ + AQ + A^ = (I + L) (D + Ay) ■ (I + U ) (D + A^)
and solve (I + L) f ■ v, (I + U ) f * v. Adding (D + A^)x = f, (D + A^)x = f
★ ★ "k
and -Aqx “ -Aqx we obtain (D + D - Ap)x = f + f - Ax = f + f -v.
J
33
Babuska's a^^orithm has some unusual round-off properties (cf. [M3])
and underflow is a problem though correctable by scaling. It handles small
row sums by taking account of dependent error distributions in the data,
and thus fits into Kahan's notion that "conserving confluence curbs ill-
condition" [Kl].
To see how the method fits in with our results, we note that, for n = 1,
since p.,q.,t i 0 we have s £ 0, s = t. + p.d'^s. , £ t.. Thus the row
J J J 1 J J J J-l J-l J
sums can only increase during the factorization. On the other hand, they
will not increase much since s.“t.+p.(s,,+q, ,)^s . £ t. + p..
J J J J-l J-l j-l J KJ
These are analogues of the results in Theorem 3.2.
3 . B. Gauss-Jordan Elimination
An algorithm closely related to block elimination is the block Gauss-
Jordan elimination
Aj = A,
Ai+i - (I + C[A.1E.)A., i = 1, . . . ,N.
This avoids back substitution in the solution of Ax = v since A^_^ is block
diagonal. The rest of the Gauss-Jordan process is
-f(D
f (1+1) - (i +
C[A. ]E.)fa)
1 ) • • •
Theorem 3.7. Let B. - B[A,]. If || B^l < 1 then ||A.+1|| £ || A. || and
II ®i+ill s II II for L =
34
Sketch of Proof. The proof parallels that of Theorem 3.2. Note = B^;
suppose || B^ || < 1. Let F^ = L[A^] + UfA^], so that A^+^ = (I - F.E.D[A ]”*)A..
= A. - F.E. + F.E.B.. This will show the first part. Now define
i ii ill.
Q = D[Ai]‘1F.E.D[A.]'1Ai, so that A.+1 = A. - D[A.]Q, D[A.+1] = D[A.](I - S) ,
S = D[Q ] . We have Q = (- B.E.)(I - B.) = - B.E. + B.E.B., so S = D[B.E.B.]
1 i l ii ill ill
and || S || £ || B.J| < 1. Also, (I - S)B.+1 =■ B. - S + Q. Let T = L - S + Q,
J = S + T, so that | J | = j S | + | T | , and Lemma 2 . 1 applies if || J j| < 1. But
J = B^(I - E^ + E^B^) , so || J 1 1 ^ || B ^ 1 1 . From Lemma 2.1 we conclude
QED
When applied to block tridiagonal matrices the block Gauss-Jordan algo-
rithm requires many more arithmetic operations than the block LU factoriza-
tion, and cannot be recommended on this basis. However, it has some numer-
ical properties that foreshadow our analysis of the more important odd-even
reduction algorithm.
The Gauss-Jordan elimination on block tridiagonals may be expressed in
the following way:
for j = 1, . . . ,N- 1
r : 1
We observe that, for 1 S i <j+k SN, -d~ 1e . , . =• u .u u , , , ,
J J,j+k J J+j. j+k-T
which is the (j,j 4- k) block of B^. Thus we are able to refine Theorem 3.7
by giving closer bounds on the first i-1 block rows of B..
l
Corollary 3.5. If ||b^|| <1 and (point) row L is contained in block row
j with 1 ^ j S i-1, then
If (point) row l is contained in block row j with i £ j £ N, then
VV = W * II Bill •
This is our first result on the convergence to zero of portions of a
matrix during the reduction to a block diagonal matrix. It also provides
a rather different situation from Gauss-Jordan elimination with partial
pivoting applied to an arbitrary matrix [PI]. In that case the multipliers
above the main diagonal can be quite large, while we have seen for block
36
tridiagonal matrices with j| B^|| < 1 that the multipliers actually decrease
in size.
3 .C . Parallel Computation
We close this section with some remarks on the applications and inher-
ent parallelism of the methods; see also [H3]. When a small number of
parallel processors are available the block and band LU factorizations
are natural approaches. This is because the solution (or reduction to tri-
angular form) of the block system d^u. = -c . is naturally expressed in
terms of vector operations on rows of the augmented matrix (d^|-Cj), and
these vector operations can be executed efficiently in parallel. However,
we often cannot make efficient use of a pipeline processor since the vec-
tor lengths might not be long enough to overcome the vector startup costs.
Our comparison of algorithms for a pipeline processor (Section 10. C) will
therefore consider the LU factorization executed in scalar mode. We also
note that the final stage of Babuska's LU-UL factorization is the solution
of a block diagonal system, which is entirely parallel.
On the other hand, the storage problems created by fill-in with these
methods when the blocks are originally sparse are well-known. For large
systems derived from elliptic equations this is often enough to disqualify
general elimination methods as competitive algorithms. The block elimina-
tion methods will be useful when the block sizes are small enough so that
they can be represented directly as full matrices, or when a stable implicit
representation can be used instead.
37
4. QUADRATIC METHODS
In Section 3 we have outlined some properties of the block LU factor-
ization and block Gauss-Jordan elimination. The properties are character-
ized by the generation of a sequence of matrices >L such that
|| B[Mi+^]|| ^ || B[1L]||. From this inequality we adopted the name "linear
methods". Now we consider some quadratic methods, which are characterized
by the inequality || B[Mi+1]|| £ || B[Mi ]2 || .
There are two basic quadratic methods, odd-even elimination and odd-
even reduction. The latter usually goes under the name cyclic reduction,
and can be regarded as a compact version of the former. We will show that
odd-even elimination is a variation of the quadratic Newton iteration for
A \ while odd-even reduction is a special case of the more general cyclic
reduction technique of Hageman and Varga [HI] and is equivalent to block
elimination on a permuted system. Both methods are ideally suited for use
on a parallel computer, as many of the quantities involved may be computed
independently of the others [H3]. Semidirect methods based on the quadratic
properties are discussed in Section 7.
4 .A. Odd-Even Elimination
We first consider odd-even elimination. Pick three consecutive block
equations from the system Ax = v, A = (a^.b^jCj):
V1V2 + bk-iVi + ck-ixk = Vi
Vk-i + Vk + Vfen “ \
ak+lXk + hk+l^+l + Ck+lXk+2 = vk+l
If we multiply equation k- 1 by
-\b
-1
k- 1 ’
equation k+1 by
■Ckbkil* and 3dd»
the
result is
38
(4.2)
("akbk-lak-l)xk-2
+ (bk ' Vk-Vk-l - Ckbk!iak+l)xk
+ (-Ckbk+lCk+l)xk+2
= (vk " - CkClVl)-
For k = 1 or N there are only two equations involved, and the necessary
modifications should be obvious. The name odd-even comes from the fact
that if k is even then the new equation has eliminated the odd unknowns
and left the even ones. This is not the only possible row operation to
achieve this effect, but others differ only by scaling and possible use of
a matrix commutativity condition.
By collecting the new equations into the block pentadiagonal system
(2)
H^x - v , it is seen that the row eliminations have preserved the fact
that the matrix has only three non-zero diagonals, but they now are farther
apart. A similar set of operations may again be applied, combining equa-
tions k-2, k and k+2 of to produce a new equation involving the unknowns
f 3)
x^_^> x^ and x^^ • These new equations form the system H^x = vv .
The transformations from one system to the next may be succinctly ex-
pressed in the following way. Let = A, v^ = v, m = flog,, N*1 , and define
Hi+1 = (I + CfV)Hi’
v(l+n - (I + C[H.])v(l), i - 1,2,. ..,m.
Block diagrams of the H sequence for N = 15 are shown in Figure 4.1a. The
distance between the three non-zero diagonals doubles at each stage until
the block diagonal matrix H , is obtained; the solution x is then obtained
nH-l
, - 1 (nH-1)
by x = H ,vv
7 nH-1
40
We note that the same elimination principle may be applied to matrices
of the periodic tridiagonal form
It is only necessary to include row operations that "wrap around" the matrix,
such as
row(N) - a b 1 row(N-l) - c b 1 row(l).
N N- 1 N I
An H sequence for N 3 16 is diagrammed in Figure 4.1b, and is constructed by
the same general rule as before, though there are serious and interesting
complications if N is not a power of 2.
Clearly, we can define an H sequence for any partitioned matrix A, and
only have to establish conditions under which the sequence may be continued.
In fact, the sequence ||b[H^]|| is quadratic when || B[A]|| < 1.
Theorem 4.1. If || B[A] || < 1 then || B[H1+1j|| * || B[H._]2 || and || H.+1 1| ^ || H . || .
If || C[A ] 1 1 x * 1 then || C[H1+1] ||x * || C^]2^ and || * || hJj .
Proof. We have 1^ - A, Hi+1 = (I + = (21 - =
2H. - H1D[H1]_1H1 - Hi(2l - D[Ht]“ 1H£) » H^I + B[H. ]) . Suppose || B[H^] || < 1,
II C[H^]||j< These assumptions are independent of one another and will
be kept separate in the proof. Since Hi * DfH^JCI - = (I - CfH^pDfH^],
we have H1+1 a D[H.](I - B[HJ2) - (I - C[H1]2)D[Hi]. Setting S - D[B[Hi]2 ] ,
41
T - D[C[Hi ]2 ] , we obtain D[Hi+1] 3 D[Hi](I - S) 3 (I - T)D[Hi] . Since
|| Sj| s || B[ HL ]2 || < 1, || T || x * || CCHi]2 || x < 1, and is invertible,
it follows that is invertible. We have B[ + ^ ] = I - D[ ] H^j3
I - (I - S)'1D[Hif1D[H.](I - B[H.]2), so B[H._]2 3 ( I - S)B[H.+1] + S;
2
similarly CfH^] ■ C[Hi+1](I - T) 4- T. Now, S is block diagonal and the
block diagonal part of B[H.+^] is null, so |B[H^] | = |(I - S)B[Hi+^]j + |s|.
2
Lemma 2.2 now applies, and we conclude || B[H^+^]|| ^ || B[H^] ||. Lemma 2.4
implies that || C[H.+1] 1^ * || CtH.]2^. To show || H.+l|| * || H.|| if
||B(H.]ff< 1, write H. + 1 3 D[H.](I- D[ H . ]" ^DfH. ] - HJBCH^) 3 D[H±J
N
+ (H. - D[H. ])B[H. ]. We then have, for 1 s 1 s n.,
i L 1 ~j=1 J
Px(H1+l) * ^(D[Ht]) + 3i(Hi - D[H.]) || B[H.]||
£ 02(D[H.]) + C£(H. - D[H.])
3 0i(H.) .
To show || H^JIj * || Hj^ if || CC^JHj < 1, write H. + 1 3 (I - C[H.1(D[H.] - H. ) x
D[H,. fSorH. ] =. D[Hi] + C[Hi](Hi - D[H^]) . We have
Y2(Hi+l) * VD[Hi]) + II C"HilH - °tHi^
5 Y£(D[H. ]) + y£(H. - D[Ht])
=* v^CH^ . QED
For a block tridiagonal matrix A we have BtH^^] 3 0 since is
block diagonal; thus the process stops at this point. The proof of Theorem
4.1 implicitly uses the convexity- like formulae
Hi+1 = Hi B[Hi] + D^Hi](I ‘
3 C[Hl Hl+ (1 - C[Ht])D[Hi].
42
As with Theorems 3.2 and 3.3, it is not difficult to find counterexamples
to show that Theorem 4.1 does not hold for arbitrary matrix norms. In
particular, it is not true that || B[A ] |J2 < 1 implies || B[H2 ] j(2 £ || B[HX] |(2
or even that || B [ H2 ] 1 1 2 < 1.
The conclusions of Theorem 4.1 are that the matrix elements, considered
as a whole rather than individually, do not increase in size, and that the
off-diagonal blocks decrease in size relative to the diagonal blocks if
j| B [A ] { j < 1. The rate of decrease is quadratic, though if || B[A]|| is very
close to 1 the actual decreases may be quite small initially. The relative
measure quantifies the diagonal dominance of H^, and shows that the diagonal
elements will always be relatively larger than the off-diagonals. One con-
clusion that is not reached is that they will be large on an absolute scale.
Indeed, consider the matrix A = (-l,2+e,-l), for which the non-zero diagonals
of are
(-l,2+et,-l) / l£:J(2+ek)
S’ ek+l = 4ek + ek '
For small e the off-diagonal elements will be about 2^ ^ and the main diagonal
2 i.
elements about 2 in magnitude. While not a serious decrease in size,
periodic rescaling of the main diagonal to 1 would perhaps be desirable.
As already hinted in the above remarks, a more refined estimate of
|| B[HX+X]|| may be made when dealing with tridiagonal matrices with constant
diagonals [S2]. Let A = (a,b,a) , so that
1
43
b-a2 b
-a2/b
b-2a2/b
-a
b-2a2/b
-a2 b
- a b
2 |a2/b I / j b-2a^/b I
-a
b-2a2.'b
0
b2-a2/b /
It follows that when 2|a| < jbj, || B[H?]|
'(2 - ||B[H1]2||), and thus |j| B^T || < |f B[H2 ] || < ||b[H1]2||. Note that
does not have constant diagonals, nor will any of the other H matrices,
but most of the rows of H, are essentially identical, and it is these rows
which determine |j B [ ] f J - This property will be important for the odd-even
reduction algorithms.
When A. = (a^,b ,Cj) is block, diagonally dominant, we can write
Hi ' (aJi),0,...,0,bj(i),0,...,0,cja)), S<D = !|(b^i))-1||(|^i)||+ |Mji> il ) ,
,(i)
,(i)
maXj j* ' Without proving it here, we can show that if 3'‘i' <1 then
3 ( i"*" 1 ) s 3^2. An irreducibility condition plus 0^ £ 1 yields the same
conclusion.
We also note that if A is symmetric then so is IL, i a 1, and thus
B[Hi]
Theorem 4.2. If A is block tridiagonal and positive definite then so is H.,
i i 1, and P(B[H.]) < 1.
Proof. First we observe that each is a 2-cyclic matrix [V3] since it can
be permuted into block tridiagonal form. Suppose that is positive definite.
Since * D[H^] is also positive definite, there exists a positive definite
matrix such that
E2t. ut so ft1+1- 2a. - a2;
44
IL and Hi+1 are positive definite iff Hi and Hi+1 are positive definite.
Since is a positive definite 2-cyclic matrix, we know from Corollary 2
of Theorem 3.6 and Corollary 1 of Theorem 4.3 of [V3] that P (B^) < 1, 5^= I - 1L
Since = I - D^IL = E^B.EL , f°lf°ws that p (B^) < 1 also. Let X be
an eigenvalue of H^. Then -1 < 1-X < 1 since PClL) < 1, and so we have
0 < X < 2 and 0<2X - X^ Si. But 2X - X^ is an eigenvalue of = 2H^-H^,
which is therefore positive definite. QED
Any quadratic convergence result such as Theorem 4.1 calls to mind the
Newton iteration. If ‘Ji(Y,Z) = Y(2I - ZY) = (2 1 - YZ)Y, Y(0) is any matrix
with || I - ZY ^||a< 1 for some norm || • ||a, and Y^i+^ = 3l(Y^ ,Z) , then
Y^ converges to Z 1 if it exists and (I - ZY^i+^) = (I - ZY^)^[H8].
It may be observed that = 3l(H^,D[H^] ^) , so while the H sequence is not
a Newton sequence there is a close relation between the two. Nor is it
entirely unexpected that Newton's method should appear, cf. Yun [Yl] and
Kung [K2] for other occurrences of Newton's method in algebraic computations.
Kung and Traub [K3] give a powerful unification of such techniques.
The Newton iteration Y^+^ * 2t(Y^ ,Z) is effective because, in a
certain sense, it takes aim on the matrix Z ^ and always proceeds towards it.
We have seen that odd-even elimination applied to a block tridiagonal matrix
eventually yields a block diagonal matrix. The general iteration
H.^ 3 S incorporates this idea by aiming the i— step towards
the block diagonal matrix DfH^]. Changing the goal of the iteration at each
step might be thought to destroy quadratic convergence, but we have seen
that it still holds in the sequence || B [ ] 1 1 -
4.B. Variants of Odd- Even Elimination
There are a number of variations of odd-even elimination, mostly dis-
cussed originally in terms of the odd-even reduction. For constant block
diagonals A = (a,b,a) with ab = ba, Hockney [H5] multiplies equation k- 1
by -a, equation k by b, equation k+1 by -a, and adds to obtain
('a2)\-2 + (b2 - 2a2)xk + ('a2)xk+2 = (bvk - avk-l ' avk+l}-
This formulation is prone to numerical instability when applied to the dis-
crete form of the Poisson equation, and a stabilization has been given by
Buneman (see [B20] and the discussion below).
With Hockney's use of Fourier analysis the block tridiagonal system
Ax * v is transformed into a collection of tridiagonal systems of the form
(- 1 , \,- 1) (sj) = (tj)> where \ ^ 2 and often X » 2. It was Hockney's obser-
vation that the reduction, which generates a sequence of Toeplitz tridiagonal
matrices (-1,X^,-1), could be stopped when l/X^ fell below the machine pre-
cision. Then the tridiagonal system was essentially diagonal and could be
solved as such without damage to the solution of the Poisson equation.
Stone [S2], in recommending a variation of cyclic reduction for the
parallel solution of general tridiagonal systems, proved quadratic conver-
gence of the ratios l/X^ and thus of the B matrices for the constant
diagonal case. Jordan [Jl], independently of our own work [H2], also extend-
ed Stone's result to general tridiagonals and the odd-even elimination algo-
rithm. The results given here apply to general block systems, and are not
restricted solely to the block tridiagonal case, though this will be the
major application.
46
We now show that quadratic convergence results hold for other formu-
lations of the odd-even elimination; we will see later how these results
carry over to odd-even reduction. In order to subsume several algorithms
into one, consider the sequences = A, v^ = v,
»i+i • Fi<2°rsi) -
v(l+1> = F . (2D[H. 1 -
1 L 1J 11 ’
where F^, G^ are nonsingular block diagonal matrices. The sequence H cor-
responds to the choices F. = I, G. = DfH.l .
i l • i
In the Hockney- Buneman algorithm for Poisson's equation with Dirichlet
boundary conditions on a rectangle we have
(a,b,a) , ab = ba, a
(i)>,
a, b
b,
F = diag (bv“')» G
ba) - b, b(i+» 1
.<» - a, a<i+1) =
DCH.f1,
b<i)2 - 2a(i)2,
-a
(i) 2
Sweet [Sll] also uses these choices to develop an odd-even reduction algo-
rithm. We again note that the matrices lose the constant diagonal prop-
erty for i 2 2, but symmetry is preserved. Each block of H. may be repre-
sented as a bivariate rational function of a and b; the scaling by en-
sures that the blocks of block rows j2i of are actually polynomials in
a and b. In fact, the diagonal blocks of these rows will be b^i+^ and the
non-zero off-diagonal blocks will be a^i+^.
Swarztrauber [S6] describes a Buneman-style algorithm for separable
elliptic equations on a rectangle. This leads to A = (a^.b^Cj) where the
blocks are all n x n, and a. = ry.I , b. = T + 3.1 , c. = y.I . In this
J J» J > J J n
case the blocks commute with each other and may be expressed as polynomials
47
in the tridiagonal matrix T. Swarztrauber' s odd-even elimination essenti-
ally uses = diag(LCM(b i- 1 > bj+2i-l^ » = Tn for j s 0 or
j 2 S + 1, G. = D[Hl ]” 1
Stone [S2] considers a similar algorithm for tri-
diagonal matrices (n = 1) but with = diagfb^^i-l bj+2i-l^’ Gi =
Note that symmetry is not preserved.
In all three cases it may be shown by induction that H
i+1
F F.H. ,
l 1 l+l
and v(i+1) = F^...F^v^+^\ These expressions signal the instability of
the basic method when F^ ^ I [B20]. However, in Section 2 .A we showed that
the B matrices are invariant under premultiplication by block diagonal
matrices, so that BfH^] * B[H^].
(i)
The Buneman modifications for stability represent the vv sequence by
sequences p^ , q^ with v^ = + q^. Defining = D[H. ],
Q. = D. - H., S. = D[Q.G.Q.] and assuming D.G.Q. = Q.G.D., we have
xi i li LXi ixiJ & i ixi xi i i’
,(» -
(i+1)
0, q
(1)
.CD
v,
p(l) + D^1(Q.p(i) + q(i))
q(l+1) = F.(Q.G.q(l) +S,p(l+1)).
Theorem 4.3. Let P = I, P. = B[H._n ] . . ] , i
o(I) ■ <VPi • V*-
Then p
(i)
(I - P^x,
Remark. Although not expressed in this form, Buzbee et al. [B20] prove spe-
cial cases of these formulae for the odd-even reduction applied to Poisson's
equation.
Proof. To begin the induction, * 0 = (I - P^)x, q*'^ = v = H^x
3 (®1^1 " Suppose p^ = (I - E\)x, q ^ * (D^P^ - Q^x. Then
■Ar
LCM a least common multiple of the matrices considered as polynomials in T.
43
p(l+1) = (I - Px)x + - P.)x + (D.P. - Ql)x)
= (I - P. + 5T 1Qi - D^Q.P. + D^D.P. - D^1Q.)x
= (I - B[H.]P.)x
= (I - Pi+1>*
, ( i+1) _ - (i+1)
= Hi+lX - W1 - Pi+l)X
(P>i+l " ^i+1 ' °i+l + °i+lPi+l)x
= <5i+lPi+l " WX‘
We note that when = DfH^] 1 and S = || B [ A ] 1 1 < 1 then ||P.J| =
and Pk(D.P. - Q.) * Pk(D . ) || P . || + Pk(Q.) * ^(D.) + ok(Q.) = Pk(H.) , so
|| D^I\ - Q^JI £ || || . It also follows that || x - || ^ || K|| || x|| and
|| q^"1 + Q^xJI £ || || || x|| . Thus p^ will be an approximation to x;
Buzbee [B17] has taken advantage of this fact in a semidirect method for
Poisson's equation.
/"'-I
4 .C. Odd- Even Reduction
4.C.I. The Basic Algorithm
The groundwork for the odd-even reduction algorithms has now been laid.
Returning to the original derivation of odd-even elimination in (4.1) and
(4.2), suppose we collect the even-numbered equations of (4.2) into a
linear system A*'2^2'* = w^, where
A< ' ' ('S2jb2 j-ia2j-l,b2 j" a2Jb2J-lc2j-l * °2 jb2 j+l*2J+l’"C2 j^j+l^J+l’iy
*«> - <*2)V, '
49
W< } " (V2j " a2jb2j-LV2j-l " C2jb2j+lV2j+l)N2 =
N2 = LN/2J.
(2)
This is the reduction step, and it is clearly a reduction since A is half
(1) (2) (2) (2) (2)
the size of A =* A. Once A x = w has been solved for x we can
compute the remaining components of x^ * x by the back substitution
X2j-1 = b2j-l(v2j-l ' a2j-lx2j-2 " C2j-lX2j)’ j = N/2'’
Since is block tridiagonal we can apply this procedure recursively
to obtain a sequence of systems A^x*'^ = , where A^ = A, w^ = v,
x(L) = (xj2i-l>Ni» - N, N.+1 = LN./2J. We write Avi; = (a^; ,c >
x^ = x . 2 l- 1 • The reduction stops with A^x^ = w^m\ m = log^ \j+l”
(note the change in notation) since A^ consists of a single block. At
this point the back substitution begins, culminating in x*"^ = x.
m rn+ 1“ i
Before proceeding, we note that if N = 2 -1 then = 2 - 1. The
odd-even reduction algorithms have usually been restricted to this particular
choice of N, but in the general case this is not necessary. However, for
A = (a,b,a)^ there are a number of important simplifications that can be
made if N = 2m - 1 [B20]. Typically these involve the special scalings men-
tioned earlier, which effectively avoid the inversion of the diagonal blocks
during the reduction.
It is also not necessary to continue the reduction recursively if
(k) (k) (k)
A x = wv can be solved more efficiently by some other method. This
leads to various polyalgorithms and the semidirect methods of Section 7.
By the way in which A^ has been constructed, it is clear that it
(4)
consists of blocks taken from H^. For example, in Figure 4.1a, A is
(*> = /-(T) u(i) „(*)■
50
the (8,8) block of . It follows immediately that || A
* H .
B[Hi] ||, and || C[A(i) ] * || C[H(l) ] || . As a corollary
of Theorem 4 . 1 we actually have
Corollary 4.1. If ||b[A]||<1 then || B[A( L+1) ] || ^ || B[A (l) ]2 || and
||A(l+l)||i If ||C[A]||1< 1 then || C[A( I+1> ) ||j s || C(A(1) ]2 ||j
anh ||A<1+1)||lS Ha'0^.
1 1.1 j _ r • ! j £ -m *■» / i _ a. A ✓ i \ _ / i I "* 1 _ II 1 1 “ 1
and A
Recalling the definition of Theorem 3.4, let X(A) = maxj(4||bj a_. || ||b_.
Theorem 4 .4 . If X(A) £ 1 then, with X^ = X(A^), we have
t(l+1) <; (X(l)/(2 - \(l)))2.
Proof. It will suffice to consider i = 1. We have
CT2j " b2ja2jb2j-lC2j-l + b2jC2jb2j+ia2j+l
(2) ,-l
3j = -a2jb2j-la2j-l’
c(2) = -c b'1 c
Cj C2j 2j+lC2j+l’
II CT2 j H ^ II b2 ja2 j II II b2 J- 1C2 j- 1 II + II b2j+ia2j+l
!Rl“a2 j)~ 1 II * (1 - V2)_1 = 2/(2 - X),
, I, (2)-l (2) || ,| (2)-l (2) ||
4Hbj aj llllbj-l Cj-lH
^ ^ lid " o>2J)“1|| ||(1 - ^2j.2)_1|| (II ^2 ja2 j
(II b2 j- Ia2 j- 1 II II b2 j-2C2 j-2 I' )
s 4(2(2 - X))2(V4)2 - ( V (2 - X))2.
b2jC2jH * 2(X/4)
b2 j- 1C2 j- 1 II }
Thus X(2) s (X(1)/(2 - X(1)))2.
51
This may be compared to earlier results for constant diagonal matrices,
(2) 2 2
which gives an equality of the form 3V = a / (2 - 3 ) . See also Theorem
4.8 of [ V3 ] .
4.C.2. Relation to Block Elimination
By regarding odd-even reduction as a compact form of odd-even elimina-
tion we are able to relate it to the Newton iteration for a"*. It is also
possible to express the reduction algorithm as block elimination on a per-
T
muted system (PAP ) (Px) = (Pv) (cf. [311], [E3], [W2]). In Section 5 we
show how this formulation also yields quadratic convergence theorems.
Specifically, let P be the permutation matrix which reorders the vector
T 0
( 1 ,2 , 3 , . . . ,N) so that the odd multiples of 2 come first, followed oy the
1 2
odd multiples of 2 , the odd multiples of 2“, etc. For N = 7,
P(l,2, . . . ,7)T = (1,3,5,7,2,6,4)T, and
This factorization gives rise to a related method, first solving L' z = Pv,
then U' (Px) 3 z. It is easily seen that z^ =
We may now apply any of the theorems about Gaussian or block elimina-
T
tion to PAP and obtain corresponding results about odd-even reduction and
the above modification. In particular, Theorem 3.2 applies though its
conclusions are weaker than those of Corollary 4.1. It follows from the
T 7 (i)
PAP = LL factorization as well as Theorem 4.2 that A will be positive
definite if A is; Theorem 4.2 also allows the conclusion that p(B[A^]) < 1.
54
The odd-even reduction is also a special case of the more general
cyclic reduction defined by Hageman and Varga [HI] for the iterative solu-
T
tion of large sparse systems. Since PAP is in the 2-cyclic form
PAP
and D? block diagonal, we can premultiply by
f 1 'A12D2
V^l0'!1
to obtain
A12D2 A21
It is clear that A
(2)
°2 ' A21°l A12>
°2 " A21°l A12 and (D1 " A12°2 A21’
0) are the rows
of ^ that have not been computed as part of the reduction. In consequence
of this formulation, if A is an M-matrix then so A
P(B[A(1+1)]) < p2(B[AU)]) < 1 [HI].
(i)
and
The block LU factorization PAP = LAU, where A is block diagonal and
L(U) is unit lower (upper) triangular, determines the reduction and back
substitution phases. Norm bounds on L and U were given in Section 3. A. 2
for general matrices with || B[A]|| < 1. For block tridiagonal matrices in
the natural ordering we have || L 1 1 ^ £ 1 + || CfAjHj , || L|| £ n|| C[A]||^ if
|| G[A] || x < 1, ||U|| S 1 + || B[A]||, || U||x * 1 + n|| B[A]|| if ||B[A]||<1,
n = maXj n ^ . For the permuted block tridiagonal system only the bounds on
|| L || and 1 1 U 1 1 ^ wi 1 1 grow, and then only logarithmically. By the construc-
tion of L and U,
55
|| L || s 1 + max || C[A(l)]|| ,
15iSm 1
INI * 1 +_ II C[A(1) ] ||,
i-1
Hull I*l+L II B[A(i> ] || ,
i-1 1
II U II s 1 + max || B[A(l) ] ||.
l^i^n
Since C[A<'l)] is block tridiagonal, p^CTA^]) 5 2n || C[A ^ ] || , so
|| C[A('I') ] || 5 2n || C[A(l) ]|| L; similarly || B[A(l) ] || 1 5 2n|| B[A(l) || . If
II B[A] || < 1, || C[A ] || < 1, we then have
|| Lllj * 1 + || CCAJHj
|| L|| 5 1 + 2nm|| C[A]||lf
|| U|| x 5 1 + 2nm|| B[A] || ,
|| U || ^1+ || N[A] || .
Much like the proof of Corollary 3.4 we also have, if || B[A]||
> x <
SO
and
lib
(i+l)
tf-
(i+l)
j
b(i)n
b2j (1 a2j )
II B[A
(i)]2
b^}|| (1 + Si> , S. = || B[A(i)]2
II (bji+1)) “1 1|
... (i+l).
cond(bj )
(b^V1!!/ (i - tj.
cond(b
(i)
2 j
)(1 + S.) / (1 - 3.).
56
Since our interest is in bounding cond(b^ji+^)
must simplify the above expression. Let = II
that
in terms of cond(b^) , we
j=1(l + 3 j ) / ( 1 - Bj), so
cond(b^1+1))
5 [max^ cond(bk)]Ki
by induction.
Using |3
we obtain
Kis"1i-i(1 + 3j)/<1- ^.i>
= (i + s.)/[(i - spn^Jd - s.)]
< (1 + 3.)/(l - S^1.
This is a rather unattractive upper bound since it may be quite large if 3^
is close to 1. Actually, the situation is not so bad, since
II II ^ Jill on(^ -i^ll _
2 j
'll<fll B[A ] || = 2^i' and
cond(b^1+1)) < [maxj^ cond(bk) ] j _ ^
2+3.
If X(A) s 1 then || <j<J} || <; X(i)/2, \(l) = X(A(l)), and
cond(b^L+1)) £ [maxk cond(bk) ] 2 + X
(i)
2 - X
(i)
^ 3 ma3C[c cond(bk>
When A is positive definite it follows from various eigenvalue- interlacing
theorems that cond(b*j^) < cond(A<'^’^) < cond(A).
4.C.3. Back Substitution
Now let us consider error propagation in the back substitution. Suppose
yki+^ ■ x,<i+1) + ?k^+^\ 1 ^ k S defines the computed solution to
. ( i+1) (i+1) „ (i+1)
A x a w
, and we compute
57
1
Then
*J- 1
(b^ )'1(W(1) - a' ' y'~ - . c
V 2j-l; V 2j-l 2i-l yi-1
.(i)
(i+1)
2J-1 yj-l
(i) (i+1).
2J-1 yj
+ e
2 j— 1 *
y(« , y(i+1>
y2j yj
§(1)
^2 j- 1
rw
' 2j- 1
.(i)
'2j-l
(b^ )"1 a^ |(i+1) _ (b^ c(i) e(i+1^
V 2j-l; a2j-l *J-1 ^j-l* C2j-1
+ e
(i)
2 j- 1 *
-(i) , -(i+1)
a2 j -j
Now define the vectors
.(i)
.(i)
zl " <ei * ° ’ «3 » °’
.(i+1)
(i+1)
(0, 5 > o, .
We have |j z^\\ - || *(i+^|| by construction, and
?(l) = B(l)z2 + zx + z2, B(l) =■ B[A(l)].
But in consequence of the construction we actually have
If || B[A] || <1 then || B(1)z,|| £ || B
§(L) || ^ max(|| B(l)z2 + zjl , || z2 || ) .
(i)
z2 II < II Z2II •
Suppose || ?(m) || £ e and || e(l)|| * e for all i =* l,...,m-l.
II ?(l) II s ™*<il B(i) || || z2|| + || z1||, || z2 || )
* II Z2II+ II ZJI
(i+Dn I, e(i)
s (m - (i+1) + 1) e + e
“ (m - i + 1) e.
Then
(by induction on i)
58
and || || , the error in the computed solution of Ax = v, is bounded by
me, m => rlog2 N + l"! . This is certainly a satisfactory growth rate for
error propagation in the back substitution. In fact, it is a somewhat
pessimistic upper bound, and we can more practically expect that errors in
the solution of x^+^ = w^+^ will be damped out quickly owing to
the quadratic decreases in || || .
4.C.4. Storage Requirements
On the matter of fill-in and storage requirements, we note that block
elimination in the natural ordering of A, as discussed in Section 3. A, re-
quires no additional block storage since there is no block fill-in. However,
it has already been seen that odd-even reduction does generate fill-in. The
amount of fill-in is the number of off-diagonal blocks of A ^ , ...,A^m\
m
which is 2(N. - 1) < 2N. Since the matrix A requires 3N - 2 storage
~i=2 L
blocks, this is not an excessive increase. Of course, we have not cons id red
the internal structure of the blocks, as this is highly dependent on the
specific system, nor have we considered temporary storage for intermediate
compu tat ions .
We do note, however, that if the factorization need not be saved, then
some economies may be achieved. For example, in the transition from A^
(2) (2) (2) (2)
to A , the first equation of A x = w is generated from the first
three equations of A^x^ = w^ . The first and third of these must be
saved for the back substitution, but there is no reason to save the second
(2) (2) (2) (2) (2) (2)
once Av xv = wv is formed. Thus the first equation of A^ xv = wv
may overwrite the second equation of A^x^ * w^ . Similar considera-
tions hold throughout the reduction, but some modifications are needed on
the CDC STAR-100 due to its restrictions on vector addressing (see Section 10. B).
59
The special scalings of odd-even reduction for systems derived from
separable elliptic PDE's allow implicit representation of the blocks of
as polynomials or rational functions of a single matrix variable.
Since the blocks of A are quite sparse this is necessary in order to limit
intermediate storage to a small multiple of the original storage. Unfor-
tunately, there does not seem to be a similar procedure that can be applied
to nonseparable PDE's. The only recourse is to store the blocks as dense
matrices, and on current computing systems this will limit them to small
sizes .
4.C.5. Special Techniques
Now suppose that A = (-p ,t + p + q , -q ) , p ,q ,t i 0, p = q = 0;
J J J J J J J J IN
we would like to apply Babuska's tricks for the LU factorization (Section
3 .A. 4) to the factorization induced by odd-even reduction. This only par-
(2)
tially succeeds. First, we can represent A
where
(2) (2) (2) (2) (2)
(-p . ,t. +p. + q . ,-q . )
j j j j j
pj
(2)
.(2)
.(2)
qj
P2 j b2j-l P2 j— 1 *
C2j + P2 j b2j-l *"2 j- 1 + q2j b2j+l C2j+1’
q2j b2j+l q2 j+1 ’
Ck + Pk + V
(2)
Note that t^^ 5 £ b2j" tbe ot^er han^» the technique of avoiding
T
the back substitution by combining the LU and UL factorizations of PAP
_T
will not work. Let A' = PAP
(I + L) (D + U) = (I + U*) (D* + L*) , (I + L)f = Pv,
JU J. ^ ^
(D + U) (Px) = f, (I + U )f 53 Pv, (D + L ) (Px) = f , so that
(D + D* - D[A' ]) (Px) =* f + f* - (L* + D(A' ] + U) (Px) .
60
The only way to have A' = L + D[A' ] + U is if no fill-in occurs in
ie
either factorization, in which case L = L[A' ], U = U[A]. In particular
fill-in occurs in odd-even reduction, so the method does not work in its
entirety. However, it will still be advantageous to do the reductions by
developing the t^ sequences. (See [SI] for some examples.)
SchrBder and Trottenberg [SI] also consider a method of "total reduc-
tion" for matrices drawn from PDE's and finite difference approximations
on a rectangular grid. Numerical stability and further developments are
discussed in [Sla], The total reduction is similar to odd-even reduction
(which is derived in [SI] by a "partial reduction"), but just different
enough to prevent the application of our techniques. Part of the problem
seems to be that total reduction is more closely related to the PDE than
to the matrix, though the algorithm is expressible in matrix form. Various
convergence results are undoubtedly obtainable, but we have not yet deter-
mined the proper setting for their analysis.
4 .D. Parallel Computation
Finally, the utility of the odd-even methods for parallel or pipeline
computation has been demonstrated in several studies [H3]. Stone [ S2 ] ,
Lambiotte and Voigt [L2], and Madsen and Rodrigue [Ml] discuss odd-even
★
reduction for tridiagonal systems, and conclude that it is the best method
to use when N is large (>256). Jordan [Jl] considers odd-even elimination,
which is shown in [Ml] to be most useful for middle ranges of N (64 < N < 256)
on a particular pipeline computer, the CDC STAR-100. For small values of
N the sequential LU factorization is best, owing to the effect of lower order
terms engendered by the pipeline mechanism.
★
The basic method given here is a block analogue of the method examined in [L2].
61
The inherent parallelism of the odd-even methods is in formation of
. . c , (i+1) ( i+1) ( i+1) (i+1) .
the equations of A x = wv or lL+^x = vv and in the back
substitution for odd-even reduction. It is readily seen from (4.2) that
H0 and v^ may be computed from and v^ by
compute LU factors of b^, (1 £ k £ N)
solve bk[^ ck vfc] - [ak ck vk], (1 * k £ N)
- bk ' Vk- 1
- Vk+l’ (1
s k £ N)
<2>
/V
*" Vk ' ak\-l
- Vk+i’ a
£ k £ N)
42)
- -Vk-1’ (3
s k £ N)
<2)
** "CkCk+l ’ (1
^ k s N - 2)
(Special end conditions are handled by storing zeros strategically or by
(2) (2)
shortening vector lengths slightly. A ' and w ' are computed by only gen-
erating new equations for even k; the back substitution is then rather
(2)
straightforward once x is computed.
The algorithm for odd-even reduction, N = 2m-l, all blocks n X n,
* - 1, is as follows: for i = 1 m- 1
compute LU factors of b^j^ , , (0 ^ k ^ N. ,)
Z k.+ 1 l+l
62
solve b$m) x|m) = wjm)
for i = m- 1 , . . . , 1
solve
soive D2k+1 x2k+1
(w(i) - a(i) x,(i+1)
' 2k+l a2k+l *k
c
(i) (1+1)
2k+l k+1
) , (0 £ k * Ni+1)
The ith reduction and corresponding back substitution uses
2 3 2 3
n + 0(n ))N^+^ “ 12 n arithmetic operations, and the total operation
2 3 2 3
count for the algorithm is (12-jN - 12m)n + 0(n N + n ) . This may be com-
2 3 2 3
pared to 4-jNn + 0(n N + n ) operations for solution of Ax = v by the LU
factorization .
★
For a crude estimate of execution tine on a pipeline computer, we
take tn + a, t « a, as the execution time for operations on vectors of
length N, and t, r « t « cr, as the execution time for scalar operations
[H3]. The total execution time for odd-even reduction would then be
232 3 2323
roughly (12— n + 0(n ))(tN + am) + 0(n t) versus (4— n N + 0(n N + n ))t
for the LU factorization done completely in scalar mode. A qualitative
comparison shows results similar to the more precise analysis given to
tridiagonal systems: for a fixed value of n, odd-even reduction will be-
come increasingly more effective than the LU factorization as N increases.
Odd-even elimination should maintain its utility for moderate values of N;
the factor tN + am is replaced by (tN + <7)m, but data manipulation is often
much simpler than in odd-even reduction. The actual crossover points be-
tween algorithms are implementation dependent and will be considered to
a limited extent in Section 10. C.
A critical observation on the execution time of odd-even reduction is
that roughly 50^ of the total time is spent in the first reduction from
*
A more detailed timing analysis is in Section 10. C and Appendix B.
J
63
A^ = A to . Special efforts to optimize the algorithm by solving
the reduced systems differently must take this fact into account. Further-
more, in timing each reduction it is seen that roughly 60$ of the high-
order term n^N comes from the matrix multiplications used to compute A^i+^
and 30$ from solving the 2n+l linear systems for a„, , etc. Relatively
4. Kr L
little effort is used to factor the b^'s or to make back substitutions
J
having saved the factorizations.
The use of a synchronous parallel computer creates some constraints
on implementation which are not so important on a sequential computer.
The first restriction in the interests of efficiency is that the blocks of
A should all be the same size. This allows simultaneous matrix operations
(factoring, etc.) to be programmed without resorting to numerous special
cases due to varying block sizes. In addition, pivoting causes ineffi-
ciencies when the pivot selections differ between the diagonal blocks.
As with the LU factorization, the assumption j| 3[A]|| < 1 only guarantees
that no block pivoting is needed. If A is strictly diagonally dominant
or positive definite, then no pivoting within the blocks will be required.
Thirdly, the choice of data structures to represent the matrices and
vectors must conform to the machine's requirements for the parallel or
vector operations. This is an important problem on the CDC Star owing
to its limitations on vector addressing and allocation [LI]. When A is
restricted to having small blocks this is not as serious as it would be
otherwise. In the case of the Illiac IV, the problem of data communica-
tion between processing elements must also be considered.
64
4 .E. Summary
We now summarize the results of this section. The odd-even elimina-
tion and reduction algorithms may be safely applied to the solution of
certain block tridiagonal linear systems, and in some cases various impor-
tant quantities involved in the solution decrease quadratical ly to zero.
The quadratic convergence may be explained by expressing the methods as
Newton-like iterations to compute a block diagonal matrix. Variations of
the basic methods, which have important applications to PDE's, also dis-
play quadratic convergence. Parallel computation is a natural format for
these methods, although practicalities place some restrictions on the
type of system that can be solved efficiently.
65
5. UNIFICATION: ITERATION AND FILL-IN
Having presented several linear and quadratic methods, we now treat
them in a unified manner as particular instances of a general nonstationary
iteration. We are also lead to an interpretation of quadratic convergence
in terms of matrix fill-in, and conclude that the odd-even methods are the
only instances of the general iteration which display quadratic convergence.
These results apply to arbitrary matrices satisfying |j B[A]||< 1 and not
only to the block tridiagonal case.
Let us first collect the formulae for the block LU factorization,
block Gauss- Jordan, and odd-even elimination.
LU: A.+ 1 = (I - L[Ai]EiD[Ai]“1)A1, Aj = A.
GJ: A.+1 = (I - (L[.\.] + UrA.])E.D[Ai]‘1)A., Aj = A.
OE: Hi+1 = (I - (L[H.] + ^H. , ^ = A.
Now consider the function F(X,Y,Z) = (2 1 - XY)Z, and observe that
Ai+1 = F(D^Ai^ + UAi]Ei, DCA.f1, A.)
Ai+1 - F(D[A.] + (L[At] + U[Ai])Ei, DfAj'1, A.)
Hi+1 = F(Hi’ D^Hi^1» Ht)
Thus it is necessary to analyze sequences based on repeated applications
of F.
For fixed nonsingular Y we can define the residual matrices
r (M) = I - MY, s(M) * I - YM, so F(X,Y,Z) = (I + r(X))Z and
r (F (X,Y,Z) ) - r(Z) - r(X) + r(X)r(Z),
s(F(X,Y,Z)) = S (Z) - s(X) + s (X) s (Z) .
66
By taking X = Z and assuming that || r(Z1) j| < 1 or || s(Zj) Jj < 1, the sta-
tionary iteration 3 FlZ^.Y.Z^) converges quadratical ly to Y This
is, of course, the Newton iteration.
The related operator G(X,Y,Z) = Z + X(I - YZ) 3 Z + Xs (Z) has pre-
viously been studied as a means of solving linear systems [H8]. When
Zl 3 X, Zi+^ 3 G(X,Y,Z^), Bauer [B8] observes that r(jL + ^) 3 r(X)r(Zi>,
s(Z1+1) 3 s(X)s(Zi), Z. = (I - rCXiSY"1 3 Y_1(I - s(X)S, Zi+R = G^.Y.L),
and in particular Z9^ = G(Z^,Y,Z^). This last formula also leads to the
Newton iteration, for G(Z,Y,Z) * F(Z,Y,Z), s(Z9i> = s(Z.)2.
In this section we study the nonstationary iteration Z ^ - A,
Z.,. = F(X.,Y.,Z.), where Y. = D[Z,] ^ and X. consists of blocks taken
i+l ill i L iJ i
from Z^, including the diagonal. The selection of blocks for X^ is al-
lowed to change from step to step, but we shall assume that the selection
sequence is predetermined and will not be altered adaptively. (The cost
of such an adaptive procedure, even on a highly parallel computer, would
be proh ib i t ive . )
For the solution of a linear system Ax = v we take (Z^,v^) = (A,v),
(Z^+^, = (21 ■ X^Y^) (Z^,v^) , with the sequence ending in or approach-
★ ★ * 1 *
ing (Z ,v ) , where Z is block diagonal. Then x = Z v . One necessary
condition for consistency is that 21 - X^Y^ = I + ri(X^) be nonsingular,
which holds if and only if p(r^(X^)> < 1. Since r^X^) and s^(X^) 3 I - Y^X^
are similar we could also assume p(s^(X^)> < 1. Using earlier notation,
ri(Xi> 3 C[X^], s^(X^) 3 B[X^]; the notation has been changed to provide
greater generality. In the context of an iteration converging to a block
diagonal matrix, the condition || B[A]|| < 1 (implying || s^X^) || ^ ||s^(Z^) || < 1)
states that A is sufficiently close to a particular block diagonal matrix,
namely D[A].
67
We first generalize earlier results on linear methods C ||s ^ (Z J j < 1 =»
|| si+i^Zi+x^ H 5 II s^(Z^) II ^ and t,len determine conditions under which we
obtain a quadratic method (|| || < 1 =» || si+1(Zt+1> || £ || si(Zi) ||2) .
Theorem 5 . 1 In the sequence Zi+1 = F^.Y^Z^ , suppose that
{CP ,P) | 1 ^ P ^ n} C Tt C f(p ,q) | 1 ^ p,q * N} = U,
X. = L E Z.E , Y. 3 D[Z. ]_1,
i. i pio i L i J J
i
si(M) = I - YJ1.
Then || s^CZ^) j| < 1 implies
II H * H II * H H -
Proof. Suppose || s^Z^ |j < 1. Define Q by Zi+^ = Z^ - DfZ^Q, so
Q = -Y.Z. + Y.X.Y.Z. = -s.(X.) + s.(X.)s.(Z.). Since D[s. (X.) ] = 0,
S = D[Q ] = D[s1(X1)st(Z1)]. Since s.(X.) = ^ . EpS .(Z.) Eq , || s.(X.) || * ||s.(Z.) ||
and || S || £ || si(Zi) ||2 < 1. Thus D[Zi+1] = DCZ±] (I - S) , Y.+1 = (I - S)"^.,
and s^CZ^^) = (I - S)s.+^(Zi+^) + s* Lemma 2.2 applies, so that ||s^(Z^ .) || < 1
implies || s.+1(Z1+1)|| s || st(Zi+1)||. Bue s.(Z.+1) - S.CZj) + Q -
si<Zl> * si«i> + si<Xi)si<Zi> ‘ WiC,*, + ^.Vi'Wi'V’
r = U - T^ so
+ o.(L E s. (Z.)E )
j TL p iv l q'
^ Dj(si(Zi)) .
Si(zi)
Thus || s.(Z.+1) || £ || s.(Z.) || < 1.
QED.
68
To see when the method is actually quadratic, we first establish that
2
at least some blocks of s^+^(Z^+^) are bounded in norm by || s^(Z^) || ■
Corollary 5.1. With the hypotheses of Theorem 5.1,
L E s. ,(Z., . ) E I £ I s,(Z.)
T. p l+l i+l q 1 1 i i
Proof. From
(5.1) (I - S) s.+1(Z. + 1) + S = s.(Z.) - s.(X.) + s.CX^s.CZ.),
S = D[si(Xt)s.(Zi)],
and E (I- S) = (I- S)E , we obtain
P P
(I - S) L E s . , . (Z. .)E + S
T. p l+l 1+1 q
\EP si(zi,E, - \ep Sl<Xl>E«
l+ +.EP WdW
But .i(Z1)Eq - s.(X.) - 31(Xi)Eq,
so (I - S) Sj Ep s1+1(Zl+1)Eq ■ Sj Ep s.tx^ sl(Z1)Eq - S. By Lemma 2.1,
Vp S1S-1(Z1+1)E,II S 'I Vp WS1(Z1)EJ S II sl<Xl)Sl<Zl)H S 1K<V
In certain cases Corollary 5.1 provides no useful information, since
we could have 2LE s. ,(Z. .)E =0. Block elimination is one example
i ^ p l+l l+l q
where this occurs. However, it is seen that if we are to have |js^+^(Z^+^) ||
2
£ || s.(Z.) || then we must have
(5.2) llEj.Ep st+1«1+1)Eq|| 5 II II2.
This is certainly true if T^ ia the null set, in which case we are generating
69
the sequence, so we shall assume that is non-empty and show that
(5.2) will not hold in general.
We set one final condition on T\ for notational purposes, though it
turns out to have practical consequences: if it is true that E Z E = 0
p l q
for any starting matrix Z^, then (p,q) must be in T.. For example, we
now write the Gauss- Jordan algorithm as
A. . = F(D[A.] + (L[A.] + U[A.])(Ij=1Ej) , D(A.]_1, A.).
i+1
Inclusion in of all blocks in Z^ known to be zero clearly will not
affect the sequence Z. , = F(X.,Y.,Z.). With the additional condition
i+l iii
on T., it follows that if (p,q) £ T. then E Z.E can be zero only bv
i i p i q j j
accident, and not by design.
Corollary 5.2. With the hypotheses of Theorem 5.1, only the choice = 6
yields a quadratic method for arbitrary starting matrices Z^.
Proof . Suppose that T1 ^ 6. Following the proof of Corollary 5.1, we
obtain
(I - s> Vp si+i<zi+i)E,
?L,E s . (Z . ) E + IL,E s.(X.)s.(Z.)E
p iv l q p i l l l q
which implies, by Lemma 2.1,
VP =i+i<zi+i)Eql
Er;Epsi(zi)Eq + s + sr;Epsi(xi,si(zi)EP'
1 1 r r
The presence of nonzero 1st order terms (2^,, EpS ^ (Z^)E^) in general dominates
the 2nd order terms (S + £_,E s . (X.)s . (Z.)E ), and we cannot conclude that
T'. p i i i l q
(5.2) will hold for arbitrary Z^ satisfying || s ^ (Z ^) | j < 1.
QED
70
Now let us consider what happens when block fill-in occurs, as it
will during the odd-even algorithms for block tridiagonals. For the off-
diagonal (p,q) block, p * q, we have E Z.E = 0 but E Z. ,E ^ 0. Since
p l q p 1+1 q
Zi+1 = 2Zi - Wi’ EpZi+lEq = EpXi('YiZi)Eq = EpXi Si<Zi)Eq’ and
II EpZi+lEqH * H EpXiH II Si(Zi)EqH- From <5a) and EpSEq = °» we have
(I - S) Ep Si+l(Zi+l)Eq = Ep Si(Xi} Si(Zi)Eq* WhiCh imPUeS II Ep Si+l(Zi+l>Eql
£ S + E s . (X.) s . (Z.)E s
11 p l l 11 q
Si(Zi)
(We cannot conclude that
II Ep Si+l(Zi+l)EqH * II Si+l(Zi+l)ll^ however->
As shown by these inequalities, fill-in produces small off-diagonal
blocks of Zi+1 when || s^(Z^) || < 1, and the corresponding blocks of
are quadratically small with respect to s^Z^). Thus one source of quadra-
tic convergence is from block fill-in during the process Z. , = F(X.,Y.,Z.).
l+l v i i l
In summary, we now have two explanations for quadratic convergence in
the odd-even methods when applied to block tridiagonal matrices. The first
is that the method is a variation of the quadratic Newton iteration for
A \ The second is that, in going from to H^_^, all the non-zero off-
diagonal blocks of are eliminated and all the non-zero off-diagonal
blocks of H. , arise from block fill-in. Each block of BTH. is there-
fore either zero or quadratic in 3[H^].
71
6. HIGHER-ORDER METHODS
For Poisson's equation, where A = (-1, b, -I), the choice N = 2™ - 1
leads to considerable simplification in the Buneman odd-even reduction
algorithm. The most important properties are A^ = (-1, b^\ -I)
(constant block diagonals) and the expression of as a polynomial in
b. However, in some applications other choices of N are much more natural.
Sweet [S10] therefore described a generalization of odd-even reduction
which could handle these cases and still preserve constant block diagonals
and a polynomial representation. Bank's generalized marching algorithm
[B5] for Poisson's equation uses one step of (very nearly the same) gen-
eralized reduction in order to control the stability of a fast marching
algorithm. This technique is analagous to multiple shooting methods for
ordinary differential equations.
In this section we analyze Sweet's method as yet another variation
of odd-even reduction, and show that it possesses higher-order convergence
properties. The algorithm actually considered here yields Sweet's method
through a suitable block diagonal scaling, and hence our results carry
through much as before. It must be emphasized, though, that the methods
of this section are not intended as practical methods for general block
tridiagonal systems. Rather, our purpose is to elaborate on the quadra-
tic properties of odd-even reduction by showing that they are special
cases of convergence properties for a family of methods.
6. A. p-Fold Reduction
The idea behind the generalization of odd-even reduction is to take
a block tridiagonal linear system involving the unknowns x^,x9,...,
and
72
to produce a new, smaller block tridiagonal system involving only the
unknowns x ,x„ . Once this system is solved, the knowns x ....
p 2p J p’
are back- substituted into the original equations to compute the remaining
unknowns. The entire process will be called p-fold reduction; odd-even
reduction is the case p = 2.
Let p S 2 be a given integer, and define = 'Jj/ pj . Starting with
(2) (2) (2)
Ax = v we generate the reduced system A x = w , where
A(2) = (aj2)> bj2)> cj2))x2 ' x<2) = <xpj>N > and the diagonal block b^2)
has dimension n.. For 1 s k s L + 1, define
PJ 2
bkp-p+l ckp-p+l
akp-p+2 ^kp-p+2 Ckp-p+2
akp- 1 kp-
When kp- 1 > N, either extend Ax = v with trivial equations or let be
redefined as a smaller matrix using the available blocks of A. Next, re-
partition A to have diagonal blocks
V V V b2p bN2p* \+i
and call the new system Ax = v. Three adjacent block rows of K are thus
(0
The p-fold reduction is obtained by performing a 2-fold reduction on A,
and the back substitution is obtained by solving the system
By taking
\p-tH-l 'kp-cH-l *kp-p+i\
v\p- 1
kp-1
i- V
"kp- 1 /
kp- p+1
I 0
kp-p+l'
c v
kp- 1 kp- 1
the reduction step simp lies to
(6.1)
■T
f’
, (2)
'k
42) = \p ‘ akp Vi ‘ "kp Vfl’
akp °kp- 1
bkp akp "kp-1 Ckp \p+l
Ckp ^kp+l
and the back substitution becomes
74
*kp-j °kp-j °kp-j *kp-p " Vkp-j ^kp’
1 * j * P-1.
Rather than computing the entire a and y vectors, a more efficient
technique is to use LU and UL processes on the A blocks:
for each k = 1,2,...,N2 + 1,
(LU process on A^)
^k, 1 = akp-p+l ’ ^k.l = bkp- p+1 ’
^,1 Vkp-pH-l’
for j = 2 .... ,p- 1
solve tAk, j- 1 Ck,j-1 Fk, j- 1 ^
'^k,j-l Ckp-p+j- 1 \,j-l^
\,j ^P-P+j Ak,j-1
"k.j bkp-fH-j akp-p+j Ck,j-1
■\,j Vkp-fH-j akp-p+j Fk,j-1
SOlve ^.p-l^p-l \p-l \p-ll
'°Tc,p-l Ckp-1 "^k.p-l"
75
for each k = 1,2,...,^,
(UL process on A^+^)
3k+l,-l " bkp+p-i; Vk+1,-1 = Ckp+p-l’
^+1,-1 Vk]>4-p-l
for j = 2, . .. ,p-l
solve ®k+l- j+l °k+l,-j+l Fk+l,-j+ll
r akp+p- j+1 vk+l j+1 ’\+l j+1 ^
B = b — c A
k+ 1 , - j kp+p-j kp+p-j H<+l,-j+l
y
Ckp+p-j Ck+l,-j+l
k+l,-j
co. — v — c F
_k+l,-j kp+p-j kp+p-j k+l,-j+l
solve ^kp+l r\pH-l^
^•^.p+l Vk+l,-p+l C°k+ 1 , - p+ 1 -
The block tridiagonal system A^x^ = w^ is now generated by (6.1) and
solved by some procedure, and we take advantage of the earlier LU process
for in the back substitution:
for each k = 1,2 N2 + 1,
^kp-l ^kp-1 \p-l *kp-p " Vkp-1 \p
for j = p-2 1
’'kp-jH-j Fk,j \,j *kp-p " Ck,j *kp-p+j+r
6 . B . Convergence Rates
As observed in Section 2, if j) B[A]||< 1 then, under the K partition-
ing for p-fold reduction, we have || B[£]|| £ || B[A]|| . By Corollary 4.1
it follows that p-fold reduction will yield || A^ || £ || ff|| = ||a|| and
II B[A( ^ ] || s || B[£]^ || s || B[A]||^. In fact, we can provide pth order
L
76
convergence, although our first result is somewhat weaker than might be
expected in comparison to earlier results.
We will actually consider a natural generalization of odd-even elim-
ination. Let Mp be a nonsingular block (2 p— l)-diagonal matrix,
M ( m , _,..., m m . , ) .
P J j'P+1 J,0’ j,p-l N
We want to pick M such that
P
M A = (r . , 0, ..., 0, s., 0, . . . ,0, t)
p J j j
C2)
MpA is a block (2p+l) -diagonal matrix; A is constructed by selecting
the intersections of the block rows p,2p,..., and block columns p,2p,...
Thus we have
(2) , (2)
a = r , b x = s
j JP j jP
c<2) = t. , W<2) - (M v). .
J JP J P JP
Now, the jth block row of M A is
P
~^k=- p+1
m. [block row (j+k) of A],
j > k
so that
rj = mj ,-P+l
sj = mj,-icj-i + mj,obj + raj,+i Vi
ti = mj,p-i cj+P-i'
We require that
mf, , c . , , + m, , b , . + m, , , a, , , - 0,
j,k-l j+k-1 j,k j+k j ,k+l j+k+1
I k | = 1 >2 , . . . ,p- 1 ,
77
where we take m. = 0, m. =0.
J.-P J.P
The special case p = 2, = (-a^b^ ^,1,-c^b^^), has been considered
in Section 4. A. For p > 2, define
d. = b.
J.-P+l J-p+1
dj,-k
‘‘j.-k
j,k
Vj,k
- b . . - a., d . ^ , ,c., ,,
j-k j-k j ,-k- 1 j-k- 1
k = p-2,...,l
aj~k d j , - k- 1 ’
k = p-2, ... ,0
= b.
j.P-1 j+P-1
— |) « p rj ra
j+k j+k j , k+ 1 j+k+ 1 ’
k = p-2 ) t * « i If
= -c d_1
j+k j,k+r
k = p-2 0 .
Then we mav cake
j,-k ^j.O “j.-l “j.-k+l
a .
1.0
m., =v. .v. ^ ... v ,
J »k j,0 j,l j , k- 1 ’
k = 1 ,2 , . . . , p- 1 .
The d , 's are derived by an LU factorization process, and the d 's
J>“K- j >+k
by a UL process. It follows from earlier results that if || C[A]||^< 1
then [j m^ s || , so the diagonals of decrease in size as one
moves away from the main diagonal.
Our first theorem partially generalizes Theorem 4.1.
78
Theorem 6.1. Let J = (-s.^r., 0, .... 0, -s.^t.), B = (-b.^a., 0, -b.^c,),
JJ J J J J J J
B =■ L + U, J = JL + Jy, II L || S a, II U || ^ a, || B || s 2m = 8, with * s
Then || J|| * 8?, ||jl|| s , || JyH * ±3>P .
Proof. Let RQ = 0, ^ = L(I - R^'V TQ = 0, Tfc = U(I - T^)'^,
k=l,...,p-l, S=R i+T , , so that
p- 1 p- 1
(I - S)JL = L (I - Rp_2)"1L(I - Rp_3)"1L ... L(I - R0)-1L,
(I - S)JIT = U(I - T „)'1U(I - T ,)-1U ... U(I - T.)-1U.
U p-Z p- j 0
2 - 1
Define '{r = 0 > t, = a (1 - i, .) , so
0 k k- 1
Kii* V llTklls V
||JLI|, II Jyll^oP/td- 2Vl)7^ (1- V],
* 0p/[2p_I(l - 2* ^Ijd- tk)].
Define 50 = °’ 51 = l’ 5k+l = 5k ' * 6k-l’ 30 6k = det(a'1>a()(k-l)x(k-l)> °*
Vk = 3 V 5k+l’ 1 ' *k = 'k+2 'k+1’ (1 " 2ip-l) “k=0(1 " V
(1 - 2 i ) 5 = 5
Define T = 5. l1 - a 5 so t = 1,
p-1 p p+1 p-1 k k+1 k-1 L
T2 ~ 1 - 2 or*" , ~k = Tk-1 - a2Tk_9. By induction, Tk(-) = 2 k+1,
d or = -2nk6k ^ £ 0. 3ut then, for m S 4, ~k(») 2 Tk(-|) , so 2P («») 2 1.
Thus || J|| s SP and || jJ| , || jJ| * j3P.
QED.
Corollary 6.1. If ||B[A]||<1, ||l[B[A]]||, || U[B[A j ] || * \ || B[A] || , then
|| B[A(2)]|| * || B[A] || P.
In generalization of Theorem 4.4, there is the very surprising
Theorem 6.2. With X(l) = X(A(l)), X(1) SI, y.(l^ = (l-./l- X(l) )/ (1+. 4- \(i)> ,
we have
79
\(i+D ^ X(i)PfJ_±-i£
Proof . It suffices to consider i = 1, so we xjill drop the superscripts.
Define e. = b ^ d. , T] . = b . ^ d , k = 0,...,p-l, so that
J,-k j-k j ,-k j > k J+*^ J ,k
= (-a, a5
j j,-l'v j-l “j-l 3,-2
),
.-1 . , ,-l -1 .
j- P+2 a j-p+2 '’j.-p+l
X (b"1 . a. ,,
J-p+1 j-p+1
(2) = r. ,
a • JP
3
t. = (-c . 71 . 1 .) (-b . 1 c
J J J,1 J+l J+l 3,2
X (b- 1 , a . ,) ,
v j+p- 1 j+p- 1
,-1 -.-1
j+P- 2 C j+p-2 j,p-l
c(2) - t. ,
J JP
h(2) - «
b , ■“ s . j
j jp
v - 1 , -1 1 , - 1
s . = b . - a. e . b. , c. , - c. T. b a
J J J J.-l J-l J-l J J.l J+l J+l
= b . (I - a.) .
J J
Now
’ 11 ^p-i11 3 1# ii 7ij!kii s (i ■ xit ^jlk+iii /A)"1- By definins ei = i.
Ek = (1 • XEk-l/4)"1> We haVe II ^Jp-kH * V Similarly, || * Efc.
This yields || c || * XE^/2, || (I - CTj)-1|| £ 2/(2 - \E ^ . To bound \(2)
observe that
4||b<2>-1a<2)ll II bP)-1c(2> "
4 I) (I - a. )
JP
j-l j-l"
'ln " \1 a. II II e:1 J| || b} , a.
jp JP 1 1 JP.-l 1 JP-1 JP
x [l ejP.-2,! ••• ” 'jJ.-pU11 11 'Wl “jp-P+11
x 11 (1 • w
-1
--1
X I 71
-1
JP-P.2 1
b” 1 c .
JP-P JP-P
JP-P»P“ 1 II II b JP- 1 C JP- 1 1
-1
jp-P.l" 11 jp-P+1 jp-P+1
s 4
( — 2-A (hp nr;E2.
V2 - XE J 4 j-l j
x p- 1
80
(2)
Thus is bounded by the above quantity. Suppose first that X - 1,
which implies that y, = 1. Then E^ = 2n ' (n+1) , the bound simplifies to
XP, and we are done with this case. Now suppose that X < 1. It may be
shown by induction that E^ 3 (1 + u.)(l - u.n) (1 - y,n+^) , yielding
(2/(2 - XE ))~ = (1 - X) ^((1 - ilP) (1 + u,P))""> and with a little effort
P” i
the bound simplifies to the desired result. QED
Note that we have not actually shown that X^ £ 1 implies ■£ 1.
This point is now to be resolved.
Corollary 6.2.
X^ < 1, and
If
(i+1)
u
s: 1 then \(i+1)
(i) p
, (i)P . ,
£ \ with strict
inequality if
Proof . The case X^ = 1 was considered in Theorem 6.2. Suppose that
0 \ < 1, again considering only the case i = 1. It is not hard to show
1 Lp/2J i
that ((1 + n),/2>P(2/(l + yP)> = (7. (X))" , where n (X) = (P.)(l ' •
p P ~i=0
Now, II (1) = 1, 7' (X) < 0 for 0 £ \ < 1, so n (X) > 1 for 0 ^ \ < 1 and
P P P
X^“^ < XP for this case. Assuming only \ £ 1, and since \(1 + y,)-/4 = u,
we have X^ S y,p(2/(l + y,P))% which implies y/^ ^ y,P . QED
6.C . Back Substitution
Error propagation in the back substitution for p-fold reduction de-
pends on the particular implementation. In the simplest form
*kp-j ^kp- j \p-j *kp-p ^kp-j Xkp ’
errors in
*kp- j
. depend on newly generated rounding errors and propagated
errors from x, and x, . We can write the computed solution as
kp-p kp
y, , = x, . + %. , to find
ykp-j Kp- j kp- j
81
'kp-j ^p-j 'kp-p 'kp-j ^kp + ekp-j"
It follows from the analysis of the Gauss- Jordan algorithm (Corollary 3.5)
and from explicit formulae for \p_j> j> that || l| £ || B[A]|f J,
II \p-jH S II B£A]|P if II BfAlll< which leads to
II Skp-jH * H »[AiirJll Vp«+ II alt’ll 5kPl,+ II V-jll-
Errors in a newly computed component therefore depend most strongly on
errors in its closest back-substituted component, and the "influence"
decreases as the distance between components increases. Moreover, we
have
II §kp-jl( s II B[A]H VpH’ II Spl^ + II V-jll*
If e is a bound on errors introduced at each stage, it follows that errors
in the final computed result will be bounded by e lo? N.
P
Slightly different results prevail in the back substitution
^p-l ^kp- 1 \p-l ^D-p ’Vp-l ^p
for j = 2 , . . . , p- 1
‘kp-j Fk,p-j Ak,p-j Xkp-p Ck,p-j ’Scp-j+l
Errors in the
computed solution
' kp-j
*kp-J + "kp-j obey the rule
5? pm
skp-l 'kp-l "kp-p " Vkp-p "kp + Skp-1
for j = 2 , . . . , p- 1
"kp-j " -Ak,p- j "kp-p " Ck,p-j ^kp-j+1 + €kp-j’
so the errors in x depend on propagation from x. and x, , , rather
KP J xp-p kp-j+1
than x^p.p and ’Stp’ This iroplies a linear growth as in the LU factorization
82
in conjunction with logarithmic growth as in odd-even reduction. If we
again suppose that || e^^jj <: e, with || ?kp_p |(, || ckp|| * ?, it follows
that
I! ?Up_ ! II * II B[£]|| ? + e,
II Skp.jH * || B[Xj || max(5, || lkp_j+1|| ) + e
which leads to || ?kp_j|| s % + (p-l)e. The end result is that errors in
the final computed result will be bounded by ep log N.
83
7. SEMIDIRECT METHODS
Briefly, a semidirect method is a direct method that is stopped be-
fore completion. In common usage, a direct method for the solution of
Ax - v computes x exactly in a finite number of steps if exact arithmetic
is used, while an iterative method computes a sequence of approximations
( i) *
x converging to x. The semidirect methods are intermediate on this
scale, producing an approximation to x in a finite number of steps,
while an additional finite number of steps could produce x exactly. During
the computation, intermediate results are generated which approximate the
solution. If the error is small enough then the algorithm is terminated
prematurely and an approximate solution is returned.
The method of conjugate gradients is an example of a semidirect
method since it can produce a good approximation long before its usual
termination [R2]. The accelerated Parallel Gauss iteration ([H4] and
Section 8.E) is another example, as is Malcolm and Palmer's convergent LU
factorization [M2] of A = (a,b,a) , a and b positive definite with o (b *a) < -7.
Here we consider semidirect methods based on the quadratic conver-
gence properties of the odd-even methods. This will generalize various
results obtained by Hockney [H5], Stone [S2] and Jordan [Jl]. Buzbee [B17]
discusses a similar technique based on the Buneman algorithm for Poisson's
equation (cf. Theorem 4.3).
The name semi- iterative has already been appropriated for other purposes
[V3], Stone [S2] initiated the name semidirect.
34
7 .A. Incomplete Elimination
Let us first consider odd-even elimination for Ax = v, which is the
process H1 = A, v(1) = v, H = (I + C[H.])H. , v(l+1)= (I + C[H.])v(l'1 .
In Theorem 4 . 1 we proved that if |j B[A]||< 1 then || B[Hi+^]|| ^ || B[H^]“ |j,
so that the off-diagonal blocks decrease quadratical ly in magnitude rela-
tive to the diagonal blocks. This motivates the approximate solution
z^^ = D[1L] ^v^\ Because we halt the process before completion, we
call this technique incomplete odd-even elimination.
Theorem 7.1.
x - z*-1'1 II / II X |
B[H.]||.
Proof . Since v^ = H.x = D[H^](I - B[H^])x, we have z^
x - B[H.]x.
D[H.]_1v(i)
QED
Thus the B matrix determines the relative error in the simple approxi-
mation z ^ ^ . If 3 = || B[A ] || < 1, 0 < e < 1 , m = ! log N” , then
(7.1) k= 1 + max(0, min(m.
log
l°g9 e
2 V log ?
))
guarantees |j B[R^]|| £ e. For example, if
a = —
f20 = 10-6, then
k = 6 and we should have N > 32 for incomplete odd-even elimination to be
2^ — 32 — o ^
applied. In fact, |j B[Vl^] || £ 3 =2 = 10 ”, so the results of the
incomplete elimination would be much better than required.
In the case n =» 1, a^ = c^ = a, b^ - b, we have || B[H^+^]|| *
II B[Hi ]2 }| / (2 - || B[H.]2|| ). Since x2/2 < x2 / (2 - x2) < x2 for 0 < x < 1,
k must be larger than
(7.2) 1 + log
log2(e/2)'
2l log2(3 2)
85
in order to have || B[H, ]|| £ e. (By induction, 3 > 2(9/2} , and if
2k_ ^
2(3/2) > e then 3 > e.) A more refined lower bound is obtained
2 2
by defining t = 3, t . ~ T (2 - t ) , and computing the largest k such
i rtf- in n
(7.3) xR > e .
Since B[1L ] is the matrix associated with the block Jacobi iteration for
H^x = v^ , when || B [ H^. ] | j is small we can use this iteration to improve
the estimate . pn fact, the choice z^ = of h ] * is the iterate
- iC
following the initial choice z ^ = 0. If z^k’^ is an approximation to
x and z(k»1-+1) = B[Hk]z(k,l) + D[Hk]’1v(k) then || x - z(k’^|| <:
|| B[H, ]^|| || x - z(k’°^ || . k and l may now be chosen to minimize computa-
2^ 1 £
tion time subject to the constraint 3 £ e. Any of the other standard
iterative methods may be used in place of the Jacobi iteration; for related
results using the theory of non-negative matrices and regular splittings,
see [HI],
We have remarked earlier that if the tridiagonal matrix A = (a., b., c.)
j j J
is irreducibly diagonally dominant and 9'
B[H^]|| then 3^ £ 1 and
^ 9<‘i'1 2 . This is a case of practical importance in differential
equations. However, the example A = (-1, 2, -1), for which 9^ = 1,
1 £ i ^ m-1, 3^ = 0, shows that strict inequality is necessary for incom-
plete elimination to be effective.
Indeed, when 3^ is very close to 1, as will be the case for most
second-order boundary value problems, the fact that we have strict in-
equality does not imply that 3V will be small for i < m. With regard
to the particular approximation z^ = D[H^] v^\ z ^ ^ depends only on
v £, j - (2i_1 - 1) £ l s j + (2 1-1 + 1) . Thus for N = 2m - 1, i < m,
86
( i.)
the center point 2^m ^ will be independent of the boundary data. Al-
though this is not the only possible approximation, complete odd-even
elimination will be necessary to ensure convergence to the continuous
solution .
7 ,B . Incomplete Reduction
The incomplete odd-even reduction algorithm is analagous to incomplete
odd-even elimination. Complete reduction generates a sequence of block
tridiagonal systems A^x^ = w , i = l,...,m, m = ~log? N+ 1~ , solves
, (m) (m) (m) , , (m) , , , „ . ,
Ax - w exactly for x , and back substitutes to finally obtain
(1) . , . , (k) (k) (k) ,
x = x . Incomplete reduction stops with Ax = w for some k be-
Ck) (k)
tween 1 and m, solves this system approximately for y *3 x , and back
substitutes y^ to obtain y = y^ rs x^ .
We remind the reader of observations made in Section 4.D, especially
that roughly half the actual work is performed in the first reduction, one
quarter in the second, and so forth. Thus the savings in incomplete reduc-
tion would be much less than those in incomplete elimination, where the
work per step is roughly constant. For more details, see Section 10. C.
In contrast to incomplete elimination, the incomplete reduction must
deal with error propagation in the back substitution. We have already
shown that roundoff errors in complete reduction are well-behaved, and
can only grow at most logarithmically. But since the approximation errors
quite possibly might be larger than the roundoff errors, a different analy-
sis is needed. We will see that a number of interesting properties result,
including no growth of errors under certain assumptions about roundoff,
*
This observation is due to an anonymous referee of [H2].
87
and a superconvergence- like effect due to damping of errors.
An initial approximation y^ may be computed from and w^ ac-
cording to our analysis of incomplete elimination. The choice
= D[A^^ ] lwOO is a simple one, and we have -y^'* [j ^ |(B[A^^]|
The parameter k should be determined both for computational economy and
for a suitable approximation. Here the quadratic decreases of B[A^ ] are
important and useful in some (but certainly not all) situations.
(k)
Regardless of the source of the approximation yv , let us suppose
y
(id-1)
(k) (k) (k) (i)
that y^ = x) + In the back substitution, y is computed from
by
y<» = y<i+1)
y2j Yj
(i) = (b(D )-l(w(i) . a(D d+1) . c(i) yCi+D) + e(i) ,
y2j-l C 2j-l; (W2j-l a2j-l yj-l C2j-1 yj ; + 2J-1*
with representing rounding errors. We then have
A(i) = (i+1)
2j
6(i) = e (i)
2j-l 2 j- 1
/h(i) x"l,,(i)
( 2j-l} a2 j- 1 >1
+c<« 5<i+1))
2j-l J
Theorem 7.2. If j| s(l) |j S (1 - || B[A(l) ] || )
'<i+1> II then |H(i) || - ||6<1+1)|t
88
As in the analysis of roundoff error propagation for complete odd-even
reduction, we have
i(i-) II /II
|| 5^|l ^ax(|| 6^' ||,
s(l)||+ || B[A(l)]|| || 5(l+1) || ).
From this expression and the hypothesis we have || 6^ || £ || |j .
But, also by construction, j| 5^ || ^ || || .
QED
This theorem requires some explanation. First suppose that no round-
ing errors are committed, so that = 0 and all the error in y^ is
.00
derived from the approximation error in y . It now follows from
|j B[ A ] || < 1 that this error will not grow during the back substitution.
From | x
00
s || X || , we conclude that ||x-y || / ||x || ^ -y ^ || ' ||
In fact, it is easy to see that while 1 1 x — y|| = || x^ - y <'k'> ||, many
vector components of x-y will be much smaller in size than J| x - y j j .
In Figure 7 . 1 we illustrate the errors by component for the tridiagonal
system (-1, 4, -l)x = v, >1 = 31, k = m-1 = 4, v chosen so that Xj = 1,
and y^ = b^k^ “w^k^ . In the absence of rounding errors, the back sub-
stitution quickly damps out approximation errors, creating a very pleasing
superconvergence- like effect.
Now suppose that rounding errors are committed and we have the upper
bound
.(i)
£ e. If
;(k)
^ ve for some moderate constant v then we
have essentially solved the system A^x^ = w ^ exactly, and the error
propagation analysis is the same as for complete reduction. That is, we
(i)
have || || ^ (k - i + v)e, 1 ^ i s k.
;(k)
The remaining case is when || || » e. Here the semidirect method
must be considered in competition with an iterative method. The hypoth-
esis of Theorem 7.2,
6^ || (1 - || B[A^ ] || ) 5 e, should be viewed as
90
an attempt to formalize the condition j| 5^ || » e . As long as the hypoth-
esis remains true, the difference (in norm) between x^ and y^ can en-
tirely be attributed to the original approximation error. As soon as the
hypothesis is violated it is possible for the errors to grow, but we will
still be able to say that || x - y|| = || 6^ || +• 0(ke) . Although the super-
convergence effect would be moderated somewhat by roundoff, its qualitative
effects will still be felt and we can expect this bound on accumulated er-
ror to be an overestimate.
7 . C. Applicability of the Methods
We close with some remarks on the practical application of the qua-
dratic semidirect methods. In Buzbee's analysis [B17] of the Buneman algo-
rithm for the differential equation (-A + d)u = f on a rectangle with
sides c ^ and c , da nonnegative constant, it is shown that if
1 2
T2 ^ 2
+
_2
] » 1
then a truncated Buneman algorithm will be useful. Geometrically this re-
quires that the rectangle be long and thin if d is small, and hence the
block tridiagonal matrix will have small blocks when properly ordered.
By directly considering the finite difference approximation to Poisson's
equation on a rectangle with equal mesh spacing in both directions it is
possible to reach a similar conclusion; this also serves as a simple model
for more general elliptic equations. With M nodes in one direction and N
in the other, M < N, the matrix A is (-1^, P^, - 1^) ^ , P = (-1, 4, -1) .
It is readily verified that || B[A]j|= 2 | j QM 1 1 , = P^1,
AD-A046 625
UNCLASSIFIED
CARNE6 I E-MELLON UNIV PITTSBURGH PA DEPT OF COMPUTER —ETC F/6 12/1
OIRECT AND ITERATIVE METHODS FOR BLOCK TRIDIAGONAL LINEAR SYSTE— ETCCU)
APR 77 D E HELLER N00014-76-C-0370
91
\k
l D + D 1
_ rn m- i.v «
-(1 ), n = 2m,
2D
j(l - n = 2nrt-l,
a
»0' l> Dl'4> Dn'4 Dn-1 “n-2 *
The following table gives j| B[A]|| for M = 1(1)6, along with estimates
on the number of odd-even reductions needed to obtain || B[A^ ]|| < 10
k such that ]| B[A^^]||< 10 ^
M
1
2
3
4
5
6
B[A ] j
2 = -5
! * -67
| = .86
If*-
If* -
If*-
lower bounds
(7.2) (7.3)
upper bound
(7.1)
8
9
10
11
These figures again lead to the conclusion that the best use of the semi-
direct methods for elliptic equations will be with long thin regions. In
the next section we consider some iterative methods which use matrix split
tings and orderings to take advantage of this property.
92
8. ITERATIVE METHODS
In this section we consider the iterative solution of a finite dif-
ference approximation to a general nonseparable elliptic equation. Of
particular interest will be methods which employ, as part of each itera-
tive step, the direct or semi-direct solution of a block tridiagonal
linear system. Continuing our emphasis on methods for a parallel computer,
we want this system to be such that the method already considered can be
used efficiently.
8 .A. Use of Elimination
Let us first note that several of our theorems about block elimination
have applications to iterative methods. Consider the system Ax = v and
the block Jacobi iteration * B[A]x^ + D[A] ^v. Suppose that
|j B [ A ] 1 1 < 1, so the iteration converges; this condition is satisfied for
many elliptic equations with suitable partitioning [V3], If the structure
of A can be simplified by an elimination step, forming A'x = v1 , then
subsequent computations will be handled more efficiently, and since
jj B[A' ]j| s jj B[A]|| the iteration will not converge any slower.
Indeed, convergence of the block Jacobi iteration will often be quite
slow. When A is also positive definite the Jacobi semi- iterative method
(J - SI or Chebychev acceleration) [V3]
x(l+1) = vi+1(B[Alx('l) 4- D[A]_1v - x(l_1))+ x(l"lV
vx - 1, v2 - 2/(2 - o ) ,
V.+1 = (1 - P2v1/4)“1, o = o(B[A]),
will converge at a much faster rate and is a natural candidate for parallel
93
computation. If we use 0 = || B[A]|| as an upper bound on p and use the
sequence of parameters
* * , ?
v1 = 1, v2 = 2/(2 - 0Z),
Vi = (1
2 * , - 1
3 v./4) \
then the J-SI iteration for A' will not be slower than that for A. By
"preprocessing" Ax = v we can expect that the actual computation time will
be smaller; unfortunately a complete proof of this claim is not at hand.
Juncosa and Mullikin [J2] discuss the effect of elimination on itera-
tive methods when A = I - M, M 2 0, m = 0, and m SI with strict
s
inequality for at least one value of r. Their particular interest for el-
liptic boundary value problems is to order the boundary nodes followed by
the interior nodes, and then to eliminate the boundary nodes. This would
be useful for domains with curved or irregular boundaries. It is shown
that, with A' = I - M' representing the system after eliminating one point,
p(M') £ p (M) . Thus the point Jacobi iteration will not be slower. Also,
if G = (I - L[M]) ^ U[M] is the matrix describing the Gauss- Seidel itera-
tion, then o(G') ^ p(G). Similar results follow from work on G-matrices
[B4]. As in the proof of Corollary 3.2 we obtain corresponding results
for block elimination, though not under the more general assumption
II B[A] || < 1.
8 . B. Splittings
The block Jacobi iteration is created by "splitting off" the block
diagonal portion of A, and using the iteration D[A]x^+^ = (D[A]-A)x^ + v.
Many other iterative methods are also defined by a splitting A = - M2 ,
MjX*'i+^ = + v. When is a block tridiagonal matrix, various direct
and semidirect methods can be used to solve the system for x^
94
A particularly effective procedure discussed by Concus and Golub
[C4] is readily adapted to parallel computation. The equation
s -7 • (a(x,y)vu) = f
on a rectangle R, a(x,y) > £) and sufficiently smooth , ^undergoes a change
~2 ~~2
of variable w(x,y) = a(x,y) u(x,y) and scaling by a to yield the equiv-
alent equation
- Aw 4- p (x , y) w = q
.11 1
~ 2 2 ~ 2
on R, p(x,y) - a A(a ), q = a f. The shifted iteration (-A + k)w ,
n-r I
= (k - p)wn + q, k a constant, is solved in the discrete form
(-A^ + KI)w^n+^ = (KI - P)v(n) + Q, where P is a diagonal matrix. Since
the system (-A^ + KI) b = c may be solved by the fast Poisson methods and
Chebychev acceleration is applicable, the overall method is highly paral-
lel and rapidly convergent.
8.C. Multi- line Orderings
There are a number of more classical iterative schemes which involve
solution of block tridiagonal systems. We begin with the multi-line itera-
tions for a rectangular grid; the ordering of unknowns and partitioning
for 2- line iteration on a 6x6 grid is shown in Figure 8.1. A diagram of
the resulting matrix is shown in Figure 8.2.
The underlying technique of the multi- line methods is to divide a
square or long and broad rectangle into a set of long and thin rectangles.
A discrete form of the elliptic equation is then solved over each subrec-
tangle as part of the iteration. By partitioning the grid into groups of
a few lines we obtain systems that can be solved well by general block
■■■■
97
methods, particularly the semidirect methods of Section 7.
The multi- line block Jacobi iteration with Chebychev acceleration
is probably the simplest method to implement on a parallel or pipeline
computer and should be considered as a benchmark for other methods. With
an nxn grid and a k-line iteration, n = kn' , the splitting A = is
such that consists of n' decoupled block tridiagonal matrices, each of
which has kxk blocks with N =* n. By storing zeros in strategic locations
may be regarded for computational purposes as one block tridiagonal
2 ,
matrix with kxk blocks and N = nn' = n ,/k. Lambiotte [LI] points out
that only log n odd-even reductions are needed because of this construc-
tion, rather than log N reductions as for the general case. When incom-
plete odd-even reduction is used, the number of lines may be chosen ac-
cording to the guidelines given in Section 7.C for Poisson's equation.
8 .D. Spiral Orderings
Another useful ordering of unknowns for parallel computation is the
single spiral order diagrammed in Figures 8.3 and 8.4. The spiral order
essentially produces a very long and thin rectangular region, which is
represented by the tridiagonal matrix forming the central portion of the
matrix in Figure 8.4. In the splitting A =■ - M0 we take to be the
2 2
tridiagonal center of A. For an nxn grid is n xn and has || B[M^]||
about thus incomplete odd-even reduction can be used quite efficiently
to solve systems involving .
We note that the single spiral order has been used with another split-
ting to yield the block peripheral method [B9], [E2]. Here the blocks are
of the periodic tridiagonal form, and the convergence rate of the periph-
eral block SOR method is shown experimentally to be similar to the 2-line
block SOR method [BIO],
98
The double spiral ordering (Figure 8.5, 8.6) is an attempt to com-
bine features of the single spiral order with advantages of a 2-cyclic
1 2
ordering. For an nXn grid, n odd, one spiral contains -(n + 2n - 1)
1 2
points and the other contains — (n - 2n + 1) points. Multi- spiral or
-peripheral orders can also be defined.
For the single and double spiral orders, when n is even the trailing
submatrix of A looks like
The *'d elements of A should be included in the implicit part of the itera-
tion, M^ . This can be done at little extra cost, and ensures that each
row and column of has at most two non-zero elements. When n is odd a
similar modification should be made.
The mapping from the natural ordering to the peripheral or single
2
spiral ordering is defined by a permutation of the vector T = (l,2,...,n )
2
into a new vector S = (S(l) ,S(2) , . . . ,S(n )). With an nxn grid, point (i,j)
is unknown number (i-l)n + j in the natural ordering. Given T, two vec-
tors forming the coordinates (i,j) of each unknown are computed by taking
quotient and remainder:
q(k) = L(k-l)/nJ, (1 £ k s n2)
i(k) = 1 + q(k), (1 £ k £ n2)
j(k) = k - nq(k) , (1 £ k £ n2) .
For simplicity let n = 2m- 1. Define
P ( i » j)
min(i,j, 2ra-i, 2m- j) , (1 £ i, j £ n) .
The point (i,j) is part of peripheral number p(i,j), 1 ^ p(i,j) ^ m.
100
Peripheral p has 8(m-p) points in it, except for peripheral m which con-
tains the single point (m,m).
The peripheral ordering is per. 1, per. 2, . .., per. m, where each
peripheral is taken in clockwise order: top, rhs , bottom, lhs. Define
P-1
f(i,j) = 8 (m-k) = 4(2m-p) (p-1) ,
‘ls-1
the number of elements contained in peripherals 1,2 , . . . ,p(i, j) - 1, and
g(i.j) = f(i,j) +
(i-p) + (j-p) +1 if i ^ j,
8 (m-p) - (i-p) - (j-p) +1 if i > j
= 4p(2m-p) + 1
-8m + (if i £ j)
- (if i > j)
(2p + i + j)
The vector S (k) = g(i,j) defines the necessary permutation, and can be
easily evaluated on a parallel or pipeline computer. The actual reordering
of a vector may be performed on the CDC STAR by using the built-in vector
permute instructions, though these should be used sparingly since they are
quite expensive. Once the vector S is generated, the elements of A and v
may be generated directly.
Of course, the implementation of an iterative method based on an
"unnatural" ordering involves much more developmental attention than is pre-
sented here. We only mean to show the potential for parallel computation.
A more complete comparison based on actual machine capabilities has not
been performed. For such a discussion of some other methods, see Lambiotte
[LI].
101
8 .E. Parallel Gauss
We close our section on iterative methods for block tridiagonal sys-
tems with some brief remarks on the Parallel Gauss iteration due to Traub
[T2] with extensions by Heller, Stevenson and Traub [H4]. For simplicity
we consider the normalized form A = (a^, I, c^) .
In the block LU factorization A = ( , I, 0) (0, d , cj) > let
D = diag(d^, . . . ,d ) . Since D = I - L[A]D Hl[A], a natural parallel itera-
tion is = I - L[A](D^ ^ ) Hj[A], where is some initial estimate
of D. This forms the first of three iterations making up the Parallel
Gauss method; the other two parts are Jacobi iterations for the L and U
block bidiagonal systems. It is shown in [H4] that, with \(A) ^ 1,
u(A) - (1 - vTrX)/(l + , we have || D - D(0 || £ u|| D - D(l_1) || .
Techniques for reducing the constant p, are given in [H4]. As Stone [S2]
points out, this linear convergence is much less attractive than the qua-
dratic convergence of odd-even reduction. However, in certain situations
(parabolic equations or the multi- line iterations discussed here) it is
necessary to solve a sequence of related problems, and this would hope-
fully provide suitable initial estimates.
Let A' = (a\, cj^ an<* suPPose that D' has already been computed.
Since d'. - d. = a, d * c. . - a'.(d'. .) ^c' , , by Theorem 3.4
J J j j-1 j-1 J J-l j“l
d'. - d . II ^ (h (• -j—
j J" v4 V1 + ,/1-V
1
2r-r) + (r~) (
l + .a-v
= 1 - jCvO-X + • ,1-V).
Thus || D ' - D || = max j || dj - d^ || 5 1 - -j(> fl- X + , ^1- X' ) , and we may take D'
as an initial approximation to D. Now, the Parallel Gauss iteration is
102
most useful when X and X' are not too close to 1, and in this case good
initial approximations will be available by the above bound on J | D ' - D | j .
When X and X' are close to 1 the bound will be large, but in this case
the iteration should probably not be used due to slow convergence.
For tridiagonal matrices we have
so that
Id1. - d
1 j ;
di
dj - dj ‘ ' Vl,dH 'j-!
- (di-irl(ij c>i ■ ai 'j-P-
1 - ./l-X
- d.
i + 777 1 j-i j-i
+
1 + - 1- X'
a c
. 1 « .v. . 1 ■
J J-l J J-l1
Let A = max . I a .c . - a’.c'. I < (X + X1 )/4 , 6 . = Id' - d | ,
J J J-l J J-l1 J 1 j j ’
a = (1 - ..T7x)/(1 + ./l-X') , b = 2/(1 + - '1- X' ) . Then *> = 0, 5. s a 5. , + b,
1 J J-l
so by induction 5^ £ (1 + a + . . . + a^-2') b < b'' (1-a) , or 2tJ (» /l-X + yl-X’ )
By the intermediate value theorem (yl- X + . ;1- X’ )/2 = .4-X’ for some \ be-
tween \ and V , so that 5^ ^ A/-/1-X . Again we have the situation where,
-'<■
if X is close to 1, a large upper bound results, but here the iteration
should not be used. When X is not so close to 1 the bound may give more
complete information than the one derived for block tridiagonal systems.
Let us now consider the utility of the Parallel Gauss iteration for
the parabolic equation ut = (au^)^, a = a(t,x) > 0, u = u(t,x), with
0 s x SI, 0 St, and u(0,x) , u(t,0) , u(t, 1) specified. Let u1 c* u(ik, jh) ,
2
k = At, h = Ax, p = k/h . The Crank-Nicolson method for advancing to the
next time step must solve the tridiagonal system
/T , P ni+l\ i+1 /T P „i. i
(I + - P ) u = (I - I p )u ,
P1
(- a
j- 1/2* °j-l/2 + aj+l/2 * * aj+l/2
/o) •
103
We therefore have
a — / i+1 , , i+1 i+1 \ /o i+1 i~\
A Daj-l/2’ 1 + P(aj-l/2 + aj+l/2) " Paj+l/2' 2);
let us now drop the superscript i+1 for convenience.
A will always be strictly diagonally dominant because a(t,x) > 0,
and we would like to know when we also have \(A) ^ 1. Using the mean value
theorem, define x. and x. by
J J
aj-i/2 + aj+i/2 - "aj-l/2 " “ °y
bj = ax((i+l)k, Xj) ,
a 3-2/2 + aj- 1/ 2 = 2a j- 1/2 ‘ h V
f? . = a ( (i+1) k, x .) .
J x J
After some algebraic manipulation it is seen that
\(A) = max.(l + (l + -rhb,)) ^
J paj- 1/ 2 2 J
x (1 + (1 - ■§ h S.))”1.
paj-l/2 2 J
In order to guarantee \(A) s 1 it is sufficient to assume th. t h and k satis-
fy ~ h|ax(t,x)| < 1. If m^ = max|ax(t,x)| this becomes the condition
m^k £ 2h.
_2
Suppose we want the more stringent condition \(A) £ (1 + c) , c £ 0.
This may be satisfied if we choose h and k such that, for each j,
2a . / + h b . ,
CPaj-l/2 S 1 + Dhbj 2>
cpa j 2 £ 1 - phH^ ' 2 .
With ra^ = max a(t,x), a sufficient condition is cpm^ s 1 • Phm^ 2, or
104
(8.1) k £ 2h2/(2cm0 + hn^) .
For c = 0, (8.1) reduces to the previous condition on h and k, while for
c = h we have k £ 2h/ + m^) = 0(h) .
Now, one of the best features of the Crank-Nicolson method is that
k may be chosen to be 0(h) . We have seen that a small value of c still
allows this choice. Conversely, if k = rh satisfies (8.1) and rm^ < 2,
then c £ h(l - rm^/2)/ (rm^) - 0(h). When we try to pick h and k so that
the Parallel Gauss method will be effective (c » h) we are necessarily
restricted to small time steps. When h and k are chosen in the usual way
(k = 0(h)), Parallel Gauss cannot be expected to be very effective since
\(A) will be close to 1. Thus the method appears to be less attractive
for this particular application than originally believed, and odd-even
reduction would be preferred. Once several reductions have been performed,
\ and u. may be sufficiently small to warrant use of Parallel Gauss and its
variations .
9. APPLICATIONS
In this section we discuss two applications of practical importance
in which the assumption || B[A]I| < 1 is not met; in the second case it is
shown how to get around this problem.
9 .A. Curve Fitting
Cubic spline interpolation of the data (x^,f^), 1 £ i £ N + 1, gives
rise to the tridiagonal system [C6, p. 238]
h.s. + 2 (h . + h . . ) s . + h . .s,,.
l i-l v l L-l 1 1-1 1+1
= 3(f[xi_1,xi]hi + f[xi>xi+1]hi_1) .
2 £ i £ N,
h. = x
l i+l
x. .
i
We have A = (h^, 2(h^ + h^ ^) , h^ and |( B[A]|| = — , which is quite satis-
factory. However, more general problems need not behave so nicely even
though they might give rise to block tridiagonal matrices.
Cox [C7] presents an algorithm for least-squares curve fitting with
piecewise polynomials which must solve the system
The v vectors contain coefficients of the Chebychev polynomial expansion
106
of the piecewise polynomial, and the X vectors contain Lagrange multipli-
ers derived through a constrained linear least-squares problem. The C/s
always have full rank, and if there are sufficiently many data points
between adjacent pairs of knots then the H/s will all be positive defi-
nite. In this case a s traightf orvard band elimination may be used well
when the number of knots and data points is large. When insufficient data
is available the H s will not have full rank and pivoting is needed in
the band elimination. It seems that the condition j| B[A]||< 1 will not
hold for any reasonable partitioning or data. Indeed, B[A] is not even .
defined for the most obvious partition.
9 .B . Finite Elements
In Section 2 .A we noted that matrices constructed from finite element
approximations to differential equations often do not satisfy assumptions
such as jj B[A]|| < 1. Thus some other technique must be used to establish
numerical stability in a direct method or convergence in an iterative
method. Usually it is shown that A is positive definite. The following
example, admittedly a simple one, shows how an easily implemented change
of variables can yield a new system A'x' = v' such that || B[A' ]|| < 1.
Even if the change of variables is not actually made, this will show that
II B[A] It*, < 1 for some norm || • ||^, and also serves to establish numerical
stability.
Let us begin with the equation -u" = 1 on the unit interval. Using
Hermite cubics, Strang and Fix [S4, p. 59] derive the 4th order accurate
finite element equations
107
, u . , , - 2u . + u . , . u' - u'
-4(-l±i r1 ti) + h - J'h - 1,
(9.1)
, u , , - u . ,
-I/_ J±i 111
5 ' 2h
2h
> + I"j - Hi ' 2ui + "i-i*
= 0
Taking the natural choice x^ = (u ^ , u^) , we obtain the block equation
aXj_1 + bxj + cXj+1 = (1,0)
The matrix A = (a, b, c) is positive definite, but || B[A]|| = — + 3 ' (4h) .
The problem is simply that the unknowns are out of scale.
A second set of unknowns, which are in scale, is
T h h T T
x. = (u. - TT u'. , u. +- u'.) aa (u. u.ll/o) . Since the Hermite cubi
j j 2 j* j 2 j J-12 j+1/2
can be derived from B-splines by coalescing adjacent nodes, these unknowns
are actually a more natural choice than the original. We now have
res
'1 -r
,1 1,
'1 0
v0 h/
u .
J
and, multiplying the first equation of (9.1) by h and the second by h,
2 T
aXj ^ + bXj + cxj+^ = (h • 1, h • 0) ,
108
The blocks are independent of h and we have |! B[A]|j = 1, with o^(B[A]) < 1
for the first and last block tows. (Actually, this assumes Dirichlet
boundary conditions ; Neumann conditions give p^(B[A]) = 1.)
For the more general equation -u" + Qu = f, Q a constant, the follow-
ing equations are derived [S4, p. 59]:
-L (.36 u.^ - 3h u'_L + 72 u . - 36 uj+1 + 3h uj+1)
+ §0(54 Uj-1 + 13 h uj.! + 312 u. + 54 u.+1 - 13h u^)
= F j *
io(3 uj-i -huj-i + 8h uj • 3 Vi -by
+ io(-13 uj-l - 3h uj-l + 8hu^ 13 uj+1 - 3 h u'+1)
= F'..
J
Again taking = (u^ “ ^ u j > Uj + ^ uj)T anc* suitably scaling the equa-
tions, they become
ax . . + bx . + cx . , = (hF.,F'.)
J-l J J+l J j
2
where now the blocks are 0(h) perturbations of the blocks in (9.2). Omit-
5 2 U
ting details, for h sufficiently small this gives || B[A]|| = 1 - — Qh + 0(h ).
109
Since rescaling and change of variables in the linear system is equiv-
alent to a change of basis in the original function subspace, the process
should be applicable to more general differential equations solved by
finite element methods. The underlying technique is to prove something
for spline approximations, and show that the property carries over in the
collapse of nodal points yielding the Hermite cubic approximations.
Strictly in terms of modifying linear systems by a block diagonal
scaling, it is natural to ask if either of the following problems are
solvable in general:
1. Given an arbitrary block tridiagonal matrix A, find a block
diagonal matrix E such that || B[AE]|| < 1.
2. In 1., if || B[A] || < 1, find E such that || B [ AE ] 1 1 < || B[A]||.
These problems are closely related to the optimal diagonal scaling of
matrices to minimize condition numbers, and as such are outside the scope
of this thesis.
110
10. IMPLEMENTATION OF ALGORITHMS
This section contains some general remarks on the choice of algorithm
for particular applications, and a more detailed comparison of algorithms
for use on a vector computer. In the comparison we assume uniform nxn
block sizes, later specialized to n = 2. All algorithms will use the same
data layout, assumed to be core-contained, and the execution time compari-
son will be in two parts: first for a very simplified and idealized model
in which the machine vectors are defined as storage locations in arithmetic
progression, and second for a model in which the machine vectors are defined
as contiguous storage locations (arithmetic progression with increment = 1) .
The first model uses information for the CDC STAR- 100 ignoring, among other
things, the fact that machine vectors on STAR must satisfy the conditions
for the second model. This double comparison allows us to pinooint the
effect of data manipulation required by the stringent requirement on vec-
tor operands. For a description of STAR facilities necessary for this dis-
cussion, see Section 2.B and Appendix A. Implementation of algorithms for
tridiagonal linear systems, n = 1, is adequately discussed elsewhere ([H4],
[Jl], [L2], [Ml], [S2]) and these references may be considered as back-
ground material for this section.
10. A. Choice of Algorithm (Generalities)
The salient point of our analysis has been that a number of algorithms
may be applied to block tridiagonal linear systems that are likely to arise
in practice. Sometimes it may be possible to take advantage of special
features in order to reduce computing costs or estimate error. When the
computing costs are high (such as one enormous problem or a large problem
solved repeatedly) an extra effort must be made to choose the most
Ill
efficient algorithm. This choice will be difficult to make in general,
but we feel that it is not at all difficult to choose good algorithms
in many cases.
Regardless of the computer involved, algorithms for three special
problems should be implemented before proceeding to algorithms for more
general problems. The first should deal with general block tridiagonal
matrices restricted to blocks of a small uniform size, say no more than
8x8. On a standard sequential computer, band or profile methods are ap-
propriate; on a k-parallel computer some form of block LU factorization
should be used, and on a vector computer a combination of odd-even reduc-
tion and LU factorization should be used. Storage allocation and access,
vector operations on a parallel or vector computer and implementation
issues in general can be greatly simplified under the block size restric-
tion. Options include implementation as a set of subroutines for specific
block sizes, insertion of tests for use as a semidirect method, and spe-
cial scalings or pivoting strategies for efficiency or stability. This
restricted class of algorithms may be applied directly to small-block
problems or used as a module in the iterative solution of large sparse
problems as indicated in Section 8. The experience gained may or may not
be useful for implementation of methods for moderate-sized blocks, which,
in our opinion, should be treated in a different manner.
The best choice of a parallel algorithm for moderate- sized blocks is
not clear at this time. These cases represent a cross-over between effi-
cient vector implementations of band/profile methods and the odd-even
methods. Either one is likely to be satisfactory. In Lambiotte's opinion
[LI] the band methods are better for the CDC STAR, but this choice is
heavily influenced by vector addressing restrictions on the STAR (cf. our
■
112
remarks in Section 10. B) . Other parallel or vector computers may not
present such problems, and a more complete analysis is needed.
The third basic class of algorithms are the fast Poisson solvers,
including generalizations. Recent studies have shown that many finite dif-
ference approximations to elliptic equations can be computed very effi-
ciently by iterative or modification methods based on the Poisson solvers.
There are several kinds of algorithms to choose from, and all have consid-
erable inherent parallelism. Other algorithms are available for elliptic
equations, but flexibility and ease of implementation suggest that Poisson-
based algorithms should be among the first to be considered. Brandt's
multi-level adaptive methods [B131 must also be seriously considered.
Finally, we note that matrix partitioning and block elimination is a
useful technique for dealing with matrices so large that secondary storage
must be used [R4]. Known as the hypermatrix scheme in other contexts, it
has been extensively exploited by the A SKA system, a structural analysis
package developed at the University of Stuttgart [F3]. Noor and Voigt [Ml]
discuss implementation of hypermatrix schemes on the CX STAR- 100. The
method works as follows: suppose a partition tt is given (see Section 2. A),
and let A^ be the matrix A partitioned accordingly. Block sizes in ASKA
are chosen to optimize transfers between storage levels, while the choice
on STAR must also try to maximize vector lengths and hence minimize start-
ud costs. The usual situation is that a few blocks of A will be in main
TT
memory, while the bulk of A^ will be in secondary memory. For example, if
tt' is another partition such that rr refines rr' , we would perform block
elimination on A^, x^, = v , . When A^ is block, tridiagonal so is A^, .
Since the data requirements for the LU or Cholesky factorization are then
highly localized, this is a natural procedure [E3]. Moreover, the diagonal
113
blocks of A^, will themselves be block tridiagonal, though possibly of low
order. As noted in Sections 2 .A and 3. A. 2, if || BfA^] || < 1 then || B[A_r,]||
£ || jfA^]!!, and pivoting will be restricted to those blocks of A , resi-
dent in main memory. This, of course, is precisely what is desired.
10. B. Storage Requirements
Data layout is an essential part of planning an algorithm' s imple-
mentation on any parallel computer, since computational speed can be
severely decreased if data cannot be accessed conveniently. On a vector
computer this means that conceptual vectors must conform to the machine's
definition of a vector operand. We consider some consequences of this
requirement for implementation on the CDC STAR- 100, assuming that the prob-
lem is core-contained.
Storage available on the STAR, for 64-bit words, consists of 256 reg-
isters and either 512K or 1M words of main memory. A sizable secondary
disc storage is also available but will not be considered here. A vector
on STAR consists of s consecutive storage locations, 1 £ s £ 64K.
Our data arrangement for block tridiagonal linear systems with uniform
nxn blocks is illustrated in fig. 10.1 for the case n = 2, N = 7. Define
m =* rlog7N" , N* = 2m, and if N < N* then extend Ax = v with N* - N trivial
equations x^ =* 0. The matrix A and vector v can then be held in three
2
storage vectors of length n N* and one of length nN*. Taking the subdiag-
2
onal blocks as an example, the n components of a, are a. . , 1 £ i s n,
k 1 j , k
lsj5n, 1 s k £ N*, and the storage vector a is defined by
a((i - l)N*n + (j - 1)N* + k) =* a , . Thus a total of (3n^ + n)N* loca-
l j , K
tions are needed to store the original system. We do not consider the
cost of setting up the storage vectors since each of our algorithms will
114
a
b
c
V
an,i = 0
bH,l
cu,i
vl,l
3 11 >2
bU,2
Cll,2
Vl,2
ail,3
bn,3
CH,3
Vl>3
*11.4
bll,4
Cll,4
Vl,4
aH,5
b.l,5
CH,5
Vl,5
all,6
bll,6
C 1 1 , 6
Vl,6
all,7
bH>7
Cll,7 = 0
Vl,7
a11.8 " 0
bll,8= 1
Cll,9 = °
Vl,8=°
a12,l = 0
bU,2
C 12 , 1
V2,l
al2,2
b12,2
C 12 ,2
Vo 9
al2 ,3
b12,3
C12,3
V9 3
al2,4
b 12 ,4
C12,4
V2,4
a!2,5
b12,5
C12,5
V2,5
a!2,6
b 12 , 6
C 12 , 6
V2,6
ai2,7
b12,7
C 12 , 7 ~ 0
V2,7
al2,S = 0
b12,8 = °
C 12 ,8 = °
a21,l =
0
b21,l
C21,l
a21,2
b21,2
C21,2
a21,3
b21,3
C21,3
a21,4
b21,4
C21,4
a21,5
b21,5
C21,5
a21,6
b21,6
C21,6
a21,7
b21,7
C21,7 =
0
a21,8
0
b21,8 = °
C21,8 =
0
a22,l =
0
b22,l
C22,l
a22,2
b22,2
C22,2
a22,3
b22,3
C22,3
a22,4
b22,4
C22,4
*22,5
b22,5
C22,5
a22,6
b22,6
C22,6
a22,7
b22,7
C22,7 =
0
a22.8 =
0
b22.9 = 1
C22,8 =
0
Fig. 10.1. Data layout for n=2, N = 7 , N* = 8 .
115
use the same arrangement and since the details can very widely between
systems. This particular data layout has been chosen because it makes the
odd-even reduction algorithm more efficient for the true STAR model, while
the other competitive algorithms (LU factorization, odd-even elimination)
are mostly unaffected by the choice of data layout, and can do without
extra storage for trivial equations.
Since the STAR forces vector lengths to be £ 64K, we must have
?
n“N* ^ 64K = 65,536. For a given value of n, this restricts the system
size as follows:
n
2
3 j
1 4
| 5
1 6
| 7
8
64R/n2
16384
7281
4096
| 2621
1826
1337
1024
max N*
16334
4096
1
4096
2048 1
1 1
1024
1
1024
1024
These will probably not be serious restrictions in most cases.
A more important restriction on system size is our assumption that the
problem and solution will remain core-contained. Assuming a 512K main
store, and supposing that temporary storage plus the program amounts to T
times the original storage, we obtain the following restrictions on N:
n
2
3
4
5
6
7
| 8
512K T- 1
18724
8738
5041
3276
2279
1702
1310
(T 4- 1) (3n + n) T = 2
12483
5825
3360
2134
1533
1134
873
T = 3
9362
4369
2520
1638
1149
851
655
max N* T ■ 1
16384
8192
4 096
2048
2048
1024
1024
T = 2
8192
4096
2048
2048
1024
1024
512
T = 3
3192
4096
2048 I
.
1024
1024
512
512
116
1
The amount of temporary storage used depends on the algorithm and the
extent to which A and v are overwritten. It can range from practically no
9
temporary storage (other than 0(n“) registers) for the LU factorization to
T sa m for an inefficient implementation of odd-even elimination. Odd-even
reduction would typically have 1^2.
We now consider the costs of data manipulation for the basic algorithms.
Major considerations are the speed of the load store unit on STAR, for reg-
ister-memory transfer of scalars, and the requirement that vectors be stored
as contiguous memory locations.
The LU factorization, executed mostly in scalar mode: The only data manip-
ulation is to load store between memory and the register file. To a con-
siderable extent this may be overlapped with arithmetic operations. A
complete analysis requires cycle counting in an assembly language program
and timings from actual runs. Our estimate of execution time (based on an
untested cycle count) is that the algorithm will probably be I/O bound.
With this in mind, we suggest four possible implementations: (Times given
are for n = 2, and represent minimal run times assuming complete overlap
between scalar arithmetic and load stores. Actual run time will certainly
be greater. For details, see Appendices A, B.)
Version 1, a one-time solution, all scalar mode
time a 50SN - 362
Version 2, factor + solve, partial vector mode
factor, version a, time a 324N + 3421 (use for N > ~ 38)
version b, time a 374N + 1769
version c, time a 420N - 230 (use for N < — 33)
solve, time a 303N + 1039
4
117
1
The three versions of factorization use increasingly fewer vector opera-
tions and increasingly more scalar operations. If more than one system
must be solved with the same matrix A, it is best to use Version 2. Our
execution time comparison in Section 10. C will use Version 1.
Odd- even elimination, executed in vector mode: An important feature of
this algorithm and any data layout like ours based on component vectors
is that the conceptual vectors all correspond to the strictest require-
ments for machine vectors, and hence no additional data manipulation is
needed. (cf. [Ml] for the tridiagonal case.)
Odd- even reduction, executed in vector mode: If the vector computer model
allows vectors to be storage locations in general arithmetic progression,
then no additional data manipulation is needed. In a true model of the
STAR, however, we need to perform two basic operations: (1) separate a
vector into its even and odd components, and (2) rejoin the even and odd
parts of a separated vector. The first is accomplished with two STAR
compress instructions, the second with a STAR merge instruction. The odd-
even separation is needed before each reduction step can be performed, and
must be applied to vectors a, b, £, and v. We illustrate the separation
2
for b in fig. 10.2. Having set up a control vector z consisting of n N*' 2
copies of the bits 10 (a single vector operation on STAR) , we execute
compress b -> bl(l) per z
b(2j - 1) - bl(j) , (lSjS |n2N*)
12
compress b -* M(— n N*) per z
b(2j) - bl(jn2N* + j), (1 * j s jn2U*).
2 112
The cost for each compress is n N* 4- -t-(— n N*) + 92 and the total cost to
o 2
118
b
z
bl
final state
1
b
h CD
11,1
11,1
11,1,,,
11,2
0
bn,3
b <X>
11,3,,,
11,3
1
bll,5
b
11,5..,
11,4
0
bll,7
b (1)
11,7 ,,,
b
1
b
b (1)
ij,k
11,5
12,1
12,1.,,
0
b
b CD i
^ k odd,
11 ,6
12,3
b12,3
•
•
•
1 S k S n*
b
b
12,7
12,3
21,1
321,2
321,3
>21,4
1
0
1
0
1
0
b
21,3
b
b
b
b
22,1
22,2
22.3
22.4
0
1
0
1
0
22,3
b
b (1)
22,5
b
°?2 s
b (1>
22,7
22,7
b ,
b (2)
11,2
ll,l.o.
->
b
b (2)
11,4
H,3.o,
b
b (2)
b (2)
11,6
b
b12jl(2)
ij ,k ’
11,8
12,3
>
k odd,
•
•
1 £ k £ N* 2
b , „ „
b ' <2)
12,8
22,3
b
b (3)
21,2
11,1.0,
b
b <3)
21,4
11,2 .o.
b2 1,6
b (3)
12,1,0,
b <3)
b
. ’ (3)
V
ij,k *
21,8
12,2
•
1 s k ^ n* ■ 4
b
b ' (3)
22,8
22,2
J
Fig. 10.2. Odd-even separation for n=2,N=7, N*=8.
119
17 2
odd-even separate a, b, c, and v is — (3n + n)N* + 736. Our initial
o
extension of A to N* block equations makes it possible to perform this
?
step using only 3 vector compress operations, independent of n“ . If A
2
had not been extended we would need either 2 (3n + n) vector operations
2
to compress the vectors one component subvector at a time, or n vector
2
operations to generate a modified control vector consisting of n copies
of z(l:N), followed by the original 3 vector compress operations.
Lambiotte [LI] rejected use of odd-even reduction for n > 1 on this basis.
The net effect of our extended data layout is to reduce the number
of vector operations needed to manipulate vectors for the STAR, and thus
to reduce vector startup costs. On the other hand, we have increased vec-
tor lengths and hence storage from N to N*, and this may be a more impor-
tant consideration in some cases. The modified control vector approach
would then be most natural despite the increased startup costs. In any
3
case, 0(n ) vector operations are needed for the arithmetic part of the
reduction step.
So far we have considered only the odd-even separation to prepare
for the first reduction. In general, the cost of the separation to pre-
17 2
pare for the ith reduction, 1 £ i £ m - 1, is — (3n + n)N* + 736,
o L
ITt = 2m+^ and the total cost for m- 1 reductions is -^-(3n“ + n) (N* - 2)
+ 736(m-l). For the case n = 2 this yields 59. 5N* + 736m - 855. We will
see shortly that this is roughly half as much time as is spent in the
arithmetic operations, so the overhead of dealing with a restrictive
machine vector definition can be considerable.
Data manipulation for the back substitution is simplified by saving
the odd- indexed equations of A^ in appropriate form for vector operations
(see fig. 10.2 for the final state of b) . Having computed x^i+^ = the
120
even components of x^^ , we compute the odd components of x^^ and merge
the two sets into x^ . The merge operation on STAR is comparatively
more expensive than the compress operation, but we are only dealing with
the solution vector and not blocks of A. The cost of the merge operation
at stage i is 3n.\\* + 123, and the total cost for m- 1 reductions is
6n(N'v - 2) + 123(m - 1) . For n = 2 this yields 12N'V + 123m - 147.
p-fold reduction: Suitable data arrangement and use of control vectors
for p-fold reduction can readily be generalized from our discussion of
odd-even reduction. The odd- even separation generalizes to separation
into p vectors, which requires p compress instructions per storage vector;
^ 2 jjyL 1 » j
the cost per reduction would be (p + — ) (3n“ + n)N.* + 363p, N.* - p‘ ,
O 1 1
m = log N 1 . The back substitution requires merging d vectors into 1,
1 1
using p-1 merge instructions at cost of 3 - — lnN.* + 123 (p- 1) to go
, (i+1) (L) _ . , , , .....
trom x to x . It is clear that these costs are minimized bv
Jacobi iterations: Data manipulation is unnecessary when A and v or BrA'
and D[A] ^v are stored according to our component subvector scheme.
10. C. Comparative Timing Analysis.
We close with an execution time comparison of algorithms for general
block tridiagonal systems with 2x2 blocks, based on timing information for
the CDC STAR- 100. This will expand the crude comparison of the LU factor-
ization and odd-even reduction given in Section 4.D, and illustrates many
earlier statements concerning the cost of various methods. Essential con-
siderations include vector startup costs, creating a substantial penalty
for operations on short vectors, and data manipulation to deal with
121
restrictive definitions of vector operands. Section 10. B began our dis-
cussion of the latter problem.
Our initial simplified analysis is based on the following conven-
tions, which allow us to consider only arithmetic costs:
1. we assume
a. vector operands may be any sequence of storage locations
in arithmetic progression,
b. pivoting is never necessary;
2 . we ignore
a. instruction overlap capabilities,
b. address calculations, loop control and termination tests,
c. load/store costs for scalar operands,
d. storage management costs, as a consequence of assumption l.a.
Our second, and more realistic analysis will demand that vector operands
consist of consecutive storage locations, so we must consider storage manage-
ment as explained in Section 10. B. Overlap, loop control and scalar load/
store will also be considered, but not address calculation or termination
of iterations. Algorithms and timing derivations are given in Appendix B;
here we present the main results and make comparisons.
Execution times for the basic direct methods are initially estimated
as follows (m * ’’log2N?):
1. the block LU factorization, using scalar operations only (LU) :
1012N - 7 84
2. odd-even elimination, using vector operations only (OEE) :
94 . 5Nm + 12999m - 93. 5N - 1901
r
122
3. odd-even reduction, using vector operations except for the final
block system (OER) : 106. 5N + 14839m - 14717.5
4. p-fold reduction, p 2 3:
d- 1
145N + 19206m (■
log2p
) - ( 1935 lp - 19579)
We can first dispose of p-fold reduction as a competitive method for gen-
eral problems since either LU or OER is always cheaper. Execution time
ratios show the benefits of using OER in preference to LU and OEE.
N = 1023
N = 8191
limit, N — 00
OER:
LU
.23
.13
.11
OER:
OEE
.24
.11
0.
The importance of vector startups is seen by comparing OER and LU,
for the latter is faster when N < 100. OER works by successively reducing
vector lengths, and at some point this becomes inefficient. A polyalgorithm
combining the two methods is straightforward: if m £ 6 then use LU, other-
wise do m-5 reductions, solve the remaining system by LU, and back substi-
tute. This reduces the time for OER to 106.5 N + 14839m - 46903.5, and
for N 3 1023 it represents a savings of 13 i.
We also note that vector startups are the dominant term for OER when
N £ 1532, and are the dominant term for OEE when N ^ 137. Of course, OEE
uses 0(N log N) operations vs. O(N) for the other methods and soon becomes
inefficient. Nevertheless, it is faster than LU for 455 £ N ^ 1023,
although it is always slower than OER.
Our conclusion is that the OER - LU polyalgorithm is the best direct
method of those considered here.
Now consider the cost of individual eliminations or reductions with
9
123
respect to the total cost of the algorithm. For odd-even elimination,
each elimination step has roughly equal cost, the last being slightly
cheaper. Thus each step costs about Ira of the total cost. For odd-even
reduction, the cost per reduction decreases as the algorithm proceeds:
cost of the kth reduction as a fraction of total cost
1 k ° 1
2
3
4
5
6
total cost
N = 1023, m = 10
00
C^l
17$
12$
3$
7$
7$
242,622
N = 3191, m 3 13
43*
22$
12$
m
4$
3$
1,050,531
(cost for kth red. = 106.5(2m-k) + 14839 when N = 2m - 1)
For small N, when startup costs are dominant, reduction costs are approxi-
mately equal, while for large N, when true arithmetic costs are dominant,
the kth reduction costs about 12 of the whole. This has immediate con-
sequences for the effectiveness of incomplete reduction, since the omitted
reductions are only a small part of the total cost when N is large.
For our comparison of semidirect and iterative methods, the goal will
be to approximate x with a relative error of l/N, or to reduce the initial
error by a factor of l/N. Thus if we perform either k odd-even elimina-
tions or k reductions, we want || B[H^+^]j| £ l/N or || B[A^k+^]|| £ I N,
assuming jj B[A]j|< 1. Theorem 7.1 is then applied to produce the desired
approximation by incomplete elimination or reduction. Taking 9 = j (B [A ] 1 1,
k k
we have || B[H^+^]||, || B[A^k+1^]j| £ » and demanding 3“ ^ l/N yields
k s log2((log2N) / (- log2 9) ) .
The execution time for k eliminations is
94.5Nk + 12999k + 10. 5N - 105(2k) + 1205,
and the time for k reductions is
124
106. 5N + 14839k - 96(2m'k) + 1287.
Typical savings for incomplete reduction over complete reduction would be:
3 = .9, pick k £ log2log2N + 2.718
N = 1023,
k
= 7,
savings ~ 12)5
N = 8191,
k
= 7,
savings = 7^
2
8 = , pick k
2 log2log2N + .774
N = 1023,
k
= 5,
savings = 25^
N = 3191,
k
= 5,
savings = 12 £
0 3 Pick k
log.
log,N
N = 1023,
k
= 4,
savings = 33--
N = 3191,
k
3 4,
savings = 16^
Comparison to the OER - LU polyalgorithm decreases the savings.
We next consider iterative methods and their combination with semi-
direct methods. One step of the block Jacobi iteration costs 22. 5N + 3031,
while one step of the block Jacobi semi- iteration costs 26. 5N + 3408 plus
overhead to estimate parameters such as p = o(B[A1). Execution time may
be reduced to 12N + 1840 (16N 4- 2217) per iteration by directly computing
B[A] and D[A] (cost = 33. 5N + 4367 = 3 iterations).
We will use the Jacobi iteration in its simplest form when only a few
steps are demanded, and Jacobi semi- iteration when many steps are needed.
The number of iterations needed by J-SI to reduce the initial error by a
factor of 1/N are estimated to be
125
The LU factorization is always faster than OEE, and is faster than
OER for N £ 312. Vector startup costs are dominant in OEE for N £ 152,
and in OER for N £ 798. The OER-LU polyalgorithm should perform m-6 reduc-
tions, solve the remaining system by LU, and back substitute. The time is
183N + 16233m - 70677, yielding a 16$ saving over OER for N = 1023 but only
a 3$ saving for N = 3191.
The execution time for odd-even elimination is allocated as follows:
!
N = 1023 N = 3191
arithmetic
96$
96$
overhead
4 $
4$
results
37 $
98$
vector startups
13$
2$
Analysis of the total time for odd-even reduction shows the greater impor-
tance of overhead operations:
I
N = 1023
arithmetic
overhead
N =* 8191
arithmetic
overhead
results
startups
33$
40$
73$
23$
/ d/
27$
56$
results
startups
52$
11$
62$
37$
1$
38$
38$
12$
126
limit, N -* =°
results
startups
/arithmetic
53 3
-
58%
(overhead
42*
-
423
1003
-
The increase in time for OER due to overhead is 37* for N = 1023, 6l3 for
N = 8191, and 723 in the limit as N -* ®. As a consequence of the data
layout of Section 10. B, overhead startup costs are quite small, but the
overhead result costs cannot be ignored and probably cannot be decreased.
The strong requirement that machine vectors consist of contiguous storage
locations leads to much "wasted" effort in execution time.
Most of the execution time of OER is actually spent in the arithmetic
and vector compress operations:
cost of individual operation as fraction of total cost
N = 1023
N = 8191
limit, N -» “
divide
12*
12*
123
multiply
43*
35*
32*
add / subtract
18*
15*
14 3
compress
203
29*
33*
merge
43
6$
73
transmit
3*
3*
3*
The compress operation assumes more
reflecting its low startup costs.
importance as
N increases, again
As noted earlier,
most of OER'
s time is spent
in the first few reduc
tions and subsequent back substitutions.
127
p = .9, 1.5m iterations,
2
p = ^, .75m iterations,
p = -j, .5m iterations.
2
For p = — , for example, J-SI is never faster than OER.
As suggested in Section 7. A, we can combine incomplete elimination or
reduction with Jacobi iterations. If k steps of elimination or reduction
are used, followed by l iterations and back substitution, k and l should
ok l
be chosen so that 3 s L N, or
(*) k + logo ^ 2 log2((log9N) / (-log9 S) ) .
The time for the OER-J method would then be
106. 5N + 14839k - 96(2m-k) + 1287 + 2(19(2m_k) + 2615),
2
and it only remains to minimize the time subject to (*) . For 3 = ■j,
N = 1023, the resulting values are k = 2, i = 5, which represent a savings
. 2
of 37/: over OER. For 3 = — , N = 8191, the resulting values are k = 3,
2=3, which represent a savings of 15$ over OER.
Of course, to determine k and i we need to know 3. While it is true
that 3 = j| B[A]j| can be computed directly from A at a cost of 46. 5N + 3881,
or from B[A] at a cost of 15N + 308, we recognize that half the rows of
B[A] are available as part of the first odd-even reduction. An estimate
of || B[A]|| can be computed from these rows at a cost of 7.5N + 308. For
large N even this cost considerably reduces any benefits gained from using
incomplete reduction, and points out the importance of an analytic estimate
of || B[A]|| rather than a computational estimate.
128
This ends our initial comparison of algorithms. We now make the
assumption that vector operands are consecutive storage locations, allow
overlap of scalar operations, and consider the cost of loop control and
scalar load store operations. The net effect is to lower the time esti-
mate for the LU factorization and raise the time estimate for odd-even
reduction. Odd-even elimination remains substantially unchanged. Rather
than simply recompute the figures presented earlier, which will be done in
a few cases, we try to indicate where the algorithms spend most of their
time.
Revised time estimates are given below. These should give a better
indication of actual runtime on STAR, but are probably all underestimates;
we note changes from the first set of timings, m = i log^N"! .
1. LU factorization, (LU) 600N - 150
(a 40 p decrease due to use of instruction overlap)
2. Odd-even elimination, (OEE) 98.5Nm 4- 13403m - 101. 5N - 2261
(a 4/» increase due to overhead costs)
3. Odd-even reduction. (OER) 183N + 16233m - 16108
(a 10-702 increase due to overhead costs; larger N incurs
larger increases).
Time ratios show the effect of decreasing the LU time and increasing the
OER time
N = 1023
N = 8191
limit, N — “
OER : LU
.54
.34
.31
OER:OEE
.32
.17
0.
129
cost of ith reduction as fraction of total cost
i =
1
2
3
4
5
6
N = 1023
33 4
19#
12 4
34
7 4
6 4
d =8191
45#
23#
12 4
6 i
2$
limit, N -* ®
5016
25 4
134
n
3 4
2 i
Of the total time, 86# is spent in the reduction phase and 14# in the back
substitution .
One of the effects of the OER-LU polyalgorithm is to reduce the impor-
tance of vector startup costs, as seen below.
N = 1023
= 3191
scalar time
14#
2 4
vector time
864
98 4
results
63*
91#
s tar tups
23 4
7 4
startup costs
are already
small not much
ilete reduction and incomplete
reduction +
results. (9
= II B[A] || ).
Incomplete reduction, savings
over OER
N =
1023
N = 3191
k
savings
k
savings
= .9 7
10#
7
5 4
f 5
i 4
CM
5
9 4
1 ‘
27#
4
134
130
Incomplete
reduction
+ Jacobi iteration.
saving over OER
N =
1023
N =
8191
k
J
savings k
l
sav in^s
3 •= .9
3
9
22 1° 4
6
8 i
2
3
2
5
36^ 2
6
1315
_1
47 ^ 2
9
1
5
4
As seen earlier, the seraidirect methods are best applied when the omitted
odd-even reductions make up a significant part of the OER time. There is
considerable advantage to using incomplete reduction plus Jacobi iterations
roughly doubling the savings obtained by simple incomplete reduction.
APPENDIX A
VECTOR COMPUTER INSTRUCTIONS
This section describes those features of the CDC STAR- 100 needed for
our comparison of algorithms in Section 10. Information given here is
derived from CDC reference manuals and timing measurements made at Lawrence
Livermore Laboratory [C2], and by no means represents a complete discussion
of the STAR'S facilities.
Execution times are for full-word (64 bit) floating point instructions,
given as the number of 40 nanosecond cycles needed for each instruction.
All times should be accurate to within 5^ except as noted, though actual
performance may be degraded somewhat due to memory conflicts.
A vector on STAR consists of N consecutive storage locations,
1 £ N £ 65,536. The vector length restriction is usually not serious, but
the need for consecutive storage is crucial. Other vector computers, the
TI ASC in particular, do not have this restriction. Both machines are
pipeline computers with the usual scalar instructions as well as vector
instructions .
Useful features of vector instructions on STAR include:
broadcast constants - Most vector instructions allow a scalar (the broad-
cast constant, held in a register) to be transmitted repeatedly to form
one or both of the vector operands. Use of this facility decreases exe-
cution time in a few cases.
control vectors - Bit vectors can modify the effect of vector instructions,
such as suppressing storage to selected components of the result vector
or determining which operand components are to be used in an operation.
132
Use in vector data manipulation (compress and merge) is explained
later.
offset - Modification to the base address of operands.
sign control - Certain vector operations allow the sign of the operands
to be modified before execution.
Individual classes of instructions follow.
scalar (register) instructions. Issue time is the number of cycles required
to initiate execution. Shortstop time is the time after issue when the
result is first available, but not yet placed in a register. The result
must be used at the precise time given; if not, the result cannot be used
until it is placed in a register. Result-to-register time is the number
of cycles after issue when the result is placed in a register. Overlapped
execution is possible, but the exact
effect
on execution
time is bes
determined by
actual testing, which we have
not done.
STAR
result
instruction
operation
issue
shortstop to
regis ter
62
add, normalized
2
6
11
66
subtract, normalized
2
6
11
6B
multiply, significant
2
10
15
6F
divide, significant
3
-
43
73
square root
3
-
71
78
register transmit
2
4
9
scalar (register) load/store Instructions. Up to three load/store instruc-
tions may be stacked in the load/store unit at any given time. The stack
is processed sequentially at a rate of 19 cycles per load and 16 cycles
133
per store minimum, assuming no memory conflicts. Once a memory bank has
been referenced it remains busy for 31 cycles; during this time further
references to the bank are delayed and the conflict will affect references
to other memory banks. Start to finish for a single load requires at
least 31 cycles, and a minimum of 66 cycles are required to store an item
and reload it.
STAR
instruction
operation
issue
result
to register
load store
unit busy
7E
load
3
28
19
7F
s tore
3
16
branch instructions . Execution time for a branch can range from 8 to 34
cycles, depending on the particular instruction, whether a branch is actu-
ally made, whether the target instruction is in the instruction stack, etc.
We will be only slightly pessimistic if we take 40 cycles to be the time
for a typical branch.
vector operation instructions. Execution time consists of an initial start-
up time followed by sequential appearance of result vector components, or
startup time followed by processing time of operands if a scalar result is
generated. Vector instructions may not be overlapped except for a negligible
part of the startup time. In the table below, N is the length of the input
vector(s) and V is the number of omitted items due to use of a
tor. Times for instructions D8 - DC may vary as much as +20^.
control vec-
134
■
STAR
vector
ins truction
operation
time
effect
82
add, normalized
+ 71
c . — a . + b .
Ill
86
subtract, normalized
jN + 71
c . — a . - b .
ill
8B
multiply, significant
N +159
c . — a . x b .
ill
SF
divide, significant
2N +167
c . «- a . + b .
ill
93
square root
2N +155
c. - *jr.
i i
98
transmit
jN + 91
c . *- a .
l i
99
absolute value
|n + 91
c . «- j a . ]
l 1 l 1
D8
maximum
6N - 2V + 95
c «- max a ,
N 1
DA
summation
6.5N - 2V + 122
c - a .
i-1 1
DB
product
6N - V + 113
c - "N a.
i=l i
N
DC
inner product
6N - V + 130
c - a b
i= 1
vector data
manipulation instructions.
The compress instruction transmi
selected components of one vector into a shorter result vector. Given a
vector a and a control vector z , the instruction essentiallv executes the
I code
j - 1
for i = 1, 2 , . . . , N
if zi = 1 then £c ^ j - j+l}.
The result vector c has length bitsum(z) . The merge instruction "shuffles"
[ two vectors into a third, according to the control vector:
135
length (c) = length(z) = length (a) + length(b)
j - l; k - 1
for i = 1, 2
if z. =
i
length (c)
1 then (V - a ; j - j+l}
else {Ci - b ; k - k+l}
In the table below, = length (a) , ^ = lengthfc), R - number of input
items broadcast from register.
STAR
instruction
BC
BD
operation
onpress a -> c per z
merge a,b - c per z
time
M1 + &2 + 92
3^ - 2R + 123
Execution times for these instructions may vary as much as +20
136
APPENDIX B
SUMMARY OF ALGORITHMS AND TIMINGS
We collect the algorithms discussed in the main text, and present
timing formulae in general terms for a vector computer. This is then
specialized to block tridiagonal systems with 2x2 blocks using timing
information for the CDC STAR- 100. A comparison of methods based on the
latter set of timings is given in Section 10. C. It must be emphasized
that our time estimates have not been verified by direct testing on STAR.
The algorithms solve Ax = v, A = (a . , b . , c.) with blocks of uni-
j j j N
form size n X n. The basic arithmetic operations are: (S = add or sub-
tract, M = multiply, D = divide)
1. factor an n x n block b
12 13 o
cost = — (n - n)D + — (2n - 3n + n) (M + S)
2. having factored b, solve bg = f, where f is n x n'
2
cost * n' (nD + (n“ - n) (M 4- S))
3. product of n x n and n x n' blocks
2 2
cost = n' (n M + (n -n)S)
4. difference of two n x n' blocks
cost = n' (nS)
Depending on context, the symbols S, M, D can mean either a single
scalar operation or a vector operation with vector length specified by
the algorithm.
For each algorithm we give a code in terms of blocks of A, since the
translation to individual scalar or vector operations is generally obvious.
137
In those cases where additional detail seems to be necessary it is pro-
vided. Two basic sets of timings are given, five formulae in all for
each algorithm:
1. counting arithemtic operations only, assuming machine vectors
can be any sequence of storage locations in arithmetic progres-
sion, and ignoring pivoting, instruction overlap, address cal-
culation, loop control, termination tests for iterations, load
store costs for scalar operands, and any other forms of storage
management
(a) for systems with n X n blocks in general terms for a
vector computer
(b) for systems with n x n blocks, with timing data for
the CDC STAR- 100
(c) for systems with 2x2 blocks, with STAR data
2. a more realistic model of STAR, requiring machine vectors to be
contiguous storage locations, allowing overlap of scalar opera-
tions, counting storage management and load'store costs, and
loop control to some extent, but still ignoring Divoting, ad-
dress calculations and termination tests for iterations
(a)
for
systems
with n
X n blocks.
with
STAR
data
(b)
for
sys terns
with 2
X 2 blocks.
with
STAR
data.
Additional notes are provided in many cases. Timings 1(c) and 2(b) are
used in Section 10. C.
A few remarks on the timings are necessary. Since we eventually want
estimates for 2x2 blocks, the only loop control that will be counted is
138
that on. the main loops. Operations on n x n matrices or n- vectors are
assumed to be expanded into straightline code.
The address calculations that have been ignored here would be exe-
cuted in scalar mode, and probably can be completely overlapped with
scalar arithmetic operations in the LU factorization. Their effect on
the odd-even methods would be to increase the 0(ml terms, m = "log^N-1 , but
this probably represents only a small relative increase since the 0(m)
terms consist of vector startup costs and are already large. While scaLar
operations can sometimes destroy the efficiency of a vector code, we be-
lieve that our basic conclusions will not be altered by including address
calculations .
139
1. Block LU factorization, using scalar operations mostly
Algorithm:
di * bi ; £i * vi
for i = 2 , ... , N
solve d. (u. g. ) = (c. f. )
l-l l-l i-l l-l l-l
d . = b . - a . u . .
l i ii-l
fi = vi - ai §1-1
solve dN xN = fN
for i = N-l, ... , 1
X. = g. - U. X.
1 1 11+1
Storage handling:
If desired we can overwrite d. -» b . , L. = a. d.V — a u -» c ,
1 l i ii-l l i l
x. -» g. - f. -* v.. The factorization loop needs to save either
1111
2
f ^ or and g^, a total of n + n stores and subsequent loads,
d., f. are saved then the second loop would be
i i r
for i = N- 1, ... , 1
solve d.x.=f.-c.x. ,
ii i i i+l
with a corresponding increase in execution time.
Variations of the basic algorithm, with estimate of minimum runt
due to load/store unit:
Version 1, a one-time solution.
initialization loads b^, v^
factorization loop loads c. ,, a., b., v.
i-l i l i
stores g^ l
initialization stores x^
d . and
i
If
ime
140
back substitution loop loads u., g.
11
stores x.
1
At 19 cycles per load and 16 per store,
2 2
total time £ {l9(4n“ 4- 2n) + 16(n“ 4- 2n) } (N- 1)
4- 19 (n2 4- n) 4- 16 (n-)
= (92n“ 4- 70n)(N-l) 4- 16n~ 4- 35n
for n = 2, total time 5 503N - 362
Version 2, factor + solve, A = ( ^ , I, 0)(0, d^ , 0 ) (0 , I, u.)
o
Factorization requires n“(3N - 2) loads for A.
Version a. factorization looo stores d. onlv, and comoutes
i
l. - a. d.\, u, = d . ^ c. by vector operations, vector length =
JJJ-ijJJ
total time £ 19n“(3N - 2) + 16n“N + -r^n" - n) D + -r(14n2 - 15n”
Z 6
D = t^N + etc.
for n = 2, total time S 323. 5N 4- 3421
Version b. factorization loop stores d^, u^, computes l. by vector
operations, vector length = N.
total time £ 19n2 (3N - 2) 4- 16n2(2N - 1)
+ \(2n~ - n)D 4- ~(8n2 - 9n~ 4- n) (M4-S)
for n = 2, total time 5 373. 5N + 1769
Version c. factorization loop stores d., u., 1..
ill
total time £ 19n2(3N - 2) + 16n2(3N - 2)
for n =* 2, total time 2 420N - 280
Solution of Ax = v, given the factorization of A.
forward elimination loads v^ initially
loop loads l., v., stores f.
i i i
solve dj gj = f j with vector operations.
N.
4- n) (H4S ) ,
141
back substitution loads initially
loop loads u., g., stores x.
i i r
2
total time s 19(2n (N-l) 4- 2nN) + 16(2nN)
+ -j(n^ + n) D + -jr(2n3 + 3n~ - 5n) (M+-S )
for n = 2, total time 2 302. 5N + 1039
Timings :
1. Version 1, for initial comparison, arithmetic operations without
overlap, load Stores ignored.
(a) t^(3n' + n)tD + ^(14n3 + 9n~ - 5n) (tM + tg ) In
- fn2tD + (2n3 + n2) (t^ + tg) }
(b) with STAR data: (70n3 + 114n~ - 2n)N - (60n3 4- 76n“)
(c) for n = 2, with STAR data: 1012N - 784
Version 1, for second comparison, allowing overlap and including
load/ stores .
In the absence of facilities for direct testing on STAR, a scalar
code was written for the special case n = 2 and time estimates
computed from that code. The details are tedious and wel 1- illustrated
by the simpler tridiagonal code in [L2], so they will not be given
here. The main time estimates are
initialize - loads + loop set up
factorization loop
Initialize
back substitution loop
total time estimated to be
200 cycles
400(N- 1)
250
200(N-1)
600N - 150 cycles
This timing is probably optimistic due to unforeseen memory and
overlap conflicts.
142
Notes :
1. Gaussian elimination, taking advantage of the block structure, using
scalar operations, requires the same execution time as the block LU
factorization .
2. If u., 1 s j sn-1, is saved in the factorization, then ||B[A ]||may
be computed following the first loop at a cost of
( (n- 1) t 4- ~ ) (nN) + (n- 1) a + c
max + max
using the vector maximum instruction. This yields 13N 4- 166 for n =
on STAR.
143
2. Odd-even elimination, using vector operations
We give the algorithm and timing for N = 2™ - 1, deriving formulae
solely in terms of N and m. Timing for other values of N is approximated
by setting m = ^og^N"! .
i- lx
vector
lengths
N, N-r
N-r
N-r
N-2r
N-2r
Algorithm: a ^ ^ = a^, etc.
for i = l,2,...,m (r = 2
for each j = 1,2,...,N
solve b.^(ry. y. c c.) = (a . ^ c.^ v.^)
J J J J J J J
take a. = y. = 0, cp. = 0 if j < 1 or j > N
J J J
(i+1) _ (i) (i) (i)
b. -b. -a. y. -c. a.,
J J J J-r j j+r
(i+l) _ (i) o (i) (i)
v. = v . -a. co . - c . cd.
J J J J-r J J+r
(i+l) _ (i)
a . = -a . or.
J J J-r
(i+l) _ (i)
c . = -c . V . ,
J J J+r
for each j = 1,2,...,N
(nH-1) _ (riH-1)
solve b . x . = v .
J J J
Storage handling: temporary vectors are needed for the factorization of
bj ^ and for , y^ , cd^ computed at each stage, and we can overwrite
. (i+l) . (i) (i+l) (i)
b. -* b . , a v -a. , etc.
J J J J
Timing :
1. arithmetic only
(a) [r(5n2 + n)T + -|(38n3 + 3n2 - 5n) t
l D b M
+ -|(38n2 - 9n2 - 5n)Tg}Nm
+ (4(5n2 + n)aD + -|(38n2 + 3n2 - 5n)a^
+ ^(38n^ - 9n“ - 5n)Cg) } m
■ {'jPn- - n)To + -|(46n3 - 3n“ + 5n)
1,3 2 i
+ ”^(46n - 27n“ + 5n)Tg } N
+ [(2n )t + (2n3 - 2n“)T + -^(n' + n) a
M S 2 D
+ -jr(-10n3 + 3n~ - 5n)<r
6 >1
+ ^r(-10n3 + 15n- - 5n)<rs]
(b) with STAR data,
£-(33n + 19n2 - n)Nm + -^(8740n3 + 2343n~ - 649n)m
-^(46n3 + n2 + n)N - |(2282n3 - 2037n2 + 649n)
(c) for n = 2, with STAR data,
94.5 Mm + 12999m - 93.5N - 1901.
2. With overwriting of a. - c^1'*, it is necessary
to develop the products one row at a time in a temporary vector, and
then use a transmit operation on the whole row. For specifics see
the discussion of odd-even reduction to follow. The added cost for
loop control and transmit instructions is
m- 1 -
40m + , = 1 2n(^n(N-21) + 91) = 40m + n"(N(m-2) 4- 1) + 182n(m-l),
or 4Nm 4- 404m - 8N - 360 for n = 2.
(a) with STAR data, total time is
£•( 38n3 + 23n2 - n)Nm + |(8740n3 4- 2343n2 4- 443n 4- 240)m
-7<46n3 + 9n2 + n)N - 7(2282n3 - 2043n2 + ’74 In)
4 6
(b) for n = 2, with STAR data, total time is
98 . 5Nm 4- 13403m W01.5N - 2261.
145
Notes :
1. || B[ A ] 1 1 may be computed in the first time through the loop at a cost
of ((2n-l)T, + T ) (nN) + ((2n-l)a, + c ). This yields 15N 4- 308
for n = 2 on STAR. || B [ H ^ ] j | may be computed at a cost of 15(N - 2m 1)
+ 303 at the appropriate time through the loop.
146
J
3. Odd-even reduction, using vector operations
As with odd-even elimination, the algorithm and timings are derived
for N = 2m - 1. For general N take n = ("log Nl in the timing formulae,
even though m = riog9N+ln is the proper choice. The algorithm is given in
detailed form, using the data layout and storage handling described in Sec-
tion 10. B. For the basic form, see Section 4.D.
Algorithm: A ^ ^ = A, = v, N* = 2°'
?
set up control vector z = (1 0 1 0 ...), length n“N* bits
cost not counted in timing, but = one instruction on STAR
for i = 1, 2, . ... m-1 (JT = 2"*’1"1' - 1, >h* = 2nH‘1'1)
,, „ ,(i) (i)
odd- even separate A ,
17 2
8 compress operations, cost = — (3n + n)N.* + 736
O 1
(i)
factor t>2j+1v , (0 j <: Ni+1
(i)
factorization overwrites b„ . ,
2j+l
need a temporary vector, length +
cost = 4(n^ - n) D + r(2n'5 - 3n- + n) (M + S) ,
i b
vector lengths = N.
soive b2j+1(l)(*2j+1(l) v2j+1
i+1
(0
®2j+l
(iV
= (a-
(i)
"2 j+1
(i)
2j+l
overwrites a
'22+1
(i)
2 j+1
<° »2j+1(1)>, (0 u s»1+1)
, etc .
need a temporary vector, length
cost 3 (2n + 1) (nD + (n2 - n) (M + S) ) ,
vector lengths = j*
- >2J<“ - -2J
<» v (1) - c
2 j- 1
2 j
(i)
(i)
2 j+1
, (i ^ j * N.+1)
overwrites b^
(i)
need a temporary vector, length
3
cost = 2n (M + S)
14 7
for simplicity, take vector lengths = N *
(i)
(i)
a„, cp.'~'-c v*' cp
O A '04—1 • TO •
(£> cp (i)
2j 92j+l *
(i+1) (i)
V =W2j -“2j +2j-l
overwrites w„ . ^
2 J
need a temporary vector, length
2
cost = 2n (M + S)
for simplicity, take vector lengths = N.
L (i+D =_*(!)* (i) ,7 < , < M .
j 2 j *2j-l ’ (2 j Ni+1)
overwrites a^^ ^ ’ Pr°dvct developed one row at a time:
for k = 1, . . . , n
for l = 1, .... n
hj ■aM.2j<i) °'n,2j-l<1)’ 5-l
for r = 2 , , n
(i)
CU ’ Cl} * *kr.2j 0r«.2].l<i)- <l s J * *1+l>
" ^t.y (1 s L s n; 1 s j SNi+i*)
3 . 3 2
cost = n M + (n -n)S, vector length »
+ n transmits, vector length a nN.+1*.
need two temporary vectors, length = nN. , and N
i+l i+1
c.(i+1) 3 -c, (° v (i), (l s j £ V . i)
j 2 j v2j+l * U J Ni+1
as for a
(i+1)
j
solve Xl(m) * w^13^ , scalar mode
cost = |(n2 + n)D + ^(2n3 + 3n2 - 5n) (M + S)
for i = m- 1, . . . , 1
x (i> = cp (i) . _ (i) _ (i+D
2 j+1 *2 j+1 *2 j+1 Xj
overwrites cd,
(i)
(i) v (i+1)
T2 j+1 xj+l
- V0<1, '"x,., , (0 * j s N. + 1)
2 j+1
cost = 2n (M + S) , vector length »
148
(i+1)
Xj -* temporary
cost = transmit, vector length = nN.+^*
merge x2j+l^’ temPorary “*
cost = 3nN.* + 123
1
Timings :
1. arithmetic only,
(a) {-j(5n2 + n)TD + -|(33n3 + 15n2 - 5n)
+ -^(38n3 + 3n2 - 5n)T }(N-1)
+ {^(5n2 + n)aQ + 'jjOSn3 + 15n2 " 5n)
+ |(33n3 + 3n2 - 5n)a }(m-l)
+ {^(n2 + n) tD + ^(2n3 + 3n2 - 5n) (t^ + tg) ]
(b) with STAR data: ^-(38n3 + 31n^ - n) (N- 1) +
4(3740n3 + 5103n2 - 649n)(m-]) +
o
(10n3 + 38n2 - 2n)
(c) for n = 2, with STAR data: 106. 5N + 14839m - 14717.5
2. with loop control and storage manipulation (compress, merge, transmit)
using STAR data,
(a) for general n, add 80m + ^-(55n2 + 43n)(N-l) + (182n + 950) (m-1)
(b) for n = 2, add 76. 5n + 1394m - 1390.5, so
total time = 183N + 16233m - 16108
L
149
Notes :
1. An estimate for || B[A]|| may be computed in the first time through the
first loop at a cost of ((2n-l)T + r ) (Nn/2) + ((2n-l)a + a ).
This yields 7.5N + 308 for n = 2 on STAR. The estimate of || B[A]|| is
simply the maximum row sum taken over every other block row. The true
value of || [A] || could be found by computing the remaining rows of B[A],
but this is more expensive and defeats the purpose of odd-even reduction.
150
4. p-fold reduction
The algorithm and timings are given for N = p™ - 1, p > 3. For general
N, take m * flog n! = m/log p.
pNi = m/log2?, m = log2N ' ; m =
log N+l is actually the
P
proper choice.
Algorithm: (see Section 6. A for details), N = d
l r
rrrt-1- i
for i = 1, .... in- 1
for each k = 1, .... Ni+1 + 1,
LU process on
for each k = 1, .... JJ. ,
l+l
UL process on ^
generate A<'i+*\
solve b^ x^ =
for i = a-1, . . . , 1
for each k = 1, .... N.. + 1
back substitute for
\p-j- 1 5 J
Timing :
2 13 2
1. (a) for general n, f(5n + n)r +T(2on + 3n - 5n) t ,
+ -j(26n3 - 3n2 - 5n)-rs3 (N-p+1)
+ [5n + n)cD + -j(26n3 + 3n2 - 5n)^
+ -|(26n3 - 3n2 - 5n) cs } (m- 1) (p- 1)
+ !^(n2 + n)tD + ^(2n2 + 3n2 - 5n)(tM+ts)}
13 2
(b) with STAR data, 2(26n + 21n - n) (N - p + 1)
+ -j(5980n + 2769n“ - 649n) (m-l)(p-l)
+ (10n3 + 38n2 - 2n)
151
(c)
using STAR data, n = 2: 145N + (19206m - 19351) (p- 1) + 228
= 145N + 19206m
(19351p - 19579)
2. data manipulation costs are outlined in Section 10. B. We have not
derived the complete timing formulae, as the method is expected to be
too slow.
Notes :
1. In the arithmetic timing, for p s 3, p-fold reduction is always cheaper
than (p-f-l)-fold reduction for any n. Comparison of 2-fold (odd-even)
and 3-fold reduction yields an N(n) such that if N s N(n)then odd-even
reduction is faster. Although we have not determined N(n) , the block
LU factorization will probably be faster than either reduction method
for N < N (n) . Similar results should hold when overhead costs are
included in the time estimates.
152
5 . Mixed odd-even reduction and LU factorization
Algorithm: do k odd-even reductions, solve A('k"rl') x('K+1^ = w<'k+^ by the
LU factorization, and back substitute.
Timing: (based on N = 2m - 1)
1 • (a) [j(5n“ + n) 4- ^(38n3 + 15n2 - 5n)
+ -|(3Sn3 + 3n2 - 5n)-r }(N + 1 - 2m-k)
1 9 13 9
+ IjCSn" + nlo^ + — (38n + 15n“ - 5n)<7^
+ |(33n3 + 3n2 - 5n);r }k
+ {— (3n + n) tQ + -|(14n3 + 9n~ - 5n)(t^ + tg)}(2m k-l)
- (n2tD + (2n3 + n2)(tM + tg)}
The timing difference between odd-even reduction and the mixed algo-
rithm will be an exoression of the form 'y (m-k) - 3 (2m_k) - v , and k
n n n
should be chosen to maximize this expression. For n = 2 using STAR
data, we have -y^ = 14839, 3^ = 905.5, v = 13023, k = m-5, yielding
(c) 106. 5N + 14839m - 46908.5
2. For n = 2 and STAR data, the polyalgorithm time is
183N + 16233k + 417(2 ‘ k) + 33, and the optimal choice is k = m-6,
yielding
(b) 183N + 16233m - 70677
15 3
6 . Incomplete reduction
Algorithm: do k odd-even reductions, solve D[A^k+^ ] y
back substitute to obtain y = y^ ^ =
(k+l) _ (k+1)
W
X.
Timing: (based on N = 2m - 1)
1. (a) {|(5n2 + n)T + -r(38n3 + 15n2 - 5n)T
^ Do M
+ x(38n3 + 3n2 - 5n)r } (N 4- 1 - 2m“k)
D S
+ {jCSn- 4- n)uD + -|(38n3 4- 15n" - 5n)
+ — (38n3 4- 3n“ - 5n)cr„}k
6 a
+ '^(n2 + n)To + |(2n3 + 3n“ - 5n) (t^ + Ts)}(2m"k - 1)
+ {^O- 4- n)a + r(2n3 4- 3n2 - 5n) (e 4- a ) }
L Do Mb
(b) general n, STAR data:
£(38n3 4- 31n2 - n)N 4- |(8740n3 4- 5103n2 - 649n) k
- (9n3 4- 6n2)(2m"k - 1) 4- -|(460n3 4- 1191n2 - 649n)
(c) using STAR data, n = 2: 106. 5N 4- 14839k - 96(2m_k) 4- 1287
2. (a) general n, STAR data:
(b)
^(38n3 4- 86n2 4- 42n)N 4- |(8740n3 4- 5103n2 4- 443n 4- 5700)k
- ^(36n3 4- 79n2 4- 51n)(2m_k - 1)
13 2
4- — (460n 4- 119 In - 1651n) 4- 80(k4-l)
n = 2, STAR data:
183N 4- 16233k - 176.5(2m'k) 4- 1113.5
and
154
7 . Jacobi iteration, serai- Iteration
basic step: = B[A]y^ + D[A]
for each j = 1, . . . , N
, , (i+1) (i) (i)
solve b.y. = v. - a.y. , - c.y. ,
J J J J J-l J J+l
Timing (identical for both vector computer models) :
1. (a) j(n2 + n) (tqN + crQ) + -|(2n3 + 15n^ - 5n)((TM + -g)N + (ct^ + cj)
1. (b) = 2. (a) using STAR data,
^(2n3 + 19n2 - n)N + |(460n3 + 3951n2 - 649n)
1. (c) = 2. (b) using STAR data, n = 2: 22. 5N 4- 3031
Optional preprocessing, times for n = 2, STAR data
(1) solve b.(a. v. cp.) =(a. c. v.),
J J J J J J J
cost = 38. 5N + 4367. The basic step is then
(i+1) _ (i) (i)
y . = co. - a.y . , - v.y . , , ,
J J J J-l J J+l
cost = 12N + 1840.
(2) factor b ^ , cost = 3.5N + 367. The basic step is unchanged from the
original, but cost is 19N + 2634. This is used in the mixed incom-
plete reduction/ Jacobi algorithm.
Semi- iteration, basic step:
y^m3 = result of Jacobi step from y^
- .al (y(m) - + y'"-1’
Vi ■ 1/(1 -
p - p(B[A]) .
p
- ^
Timing, n = 2, using STAR data: 26. 5N + 3408, or 16N + 2217 with the first
preprocessing .
Notes :
]
1. Varga [V3, p. 139] shows that for Jacobi semi- iteration ,
A conservative estimate for the number of iterations needed to obtain
|| e^ || ^ 1/n || e^ || is n = 21ogN (- log (x^- 1) ) • It follows that
the execution time for an iterative solution is 0(N log N) vs. 0(N)
for the better direct methods.
8 . Incomplete reduction plus Jacobi iterations
Algorithm: do k odd-even reductions, solve D[A^k+^ ]y(k+^ = w^k+^ f do
l Jacobi iterations to improve y(k+l)^ and back substitute.
4 , r , • , (k+1) n (k+1) (k+1) , f ‘ ,
As a result of solving D[A Jy = w , we have factors for
D[A^k+^]. This means the cost is (using STAR data, n = 2)
rn m If m. 1/
1. (c) 106. 5N + 14339k - 96(2 K) + 1287 + 2(19(2 ) + 2615).
2. (b^ 133N + 16233k - 176.5(2m"k) + 1113.5
+ 2(19(2m'k) + 2615) .
157
REFERENCES
[ B 1 ] I. Babuska, "Numerical stability in mathematical analysis," IF IP
Congress 1968, Nor th-Hol land , Amsterdam, pp. 11-23.
[B2] I. Babuska, Numerical stability in the solution of tridiagonal
matrices, Rept. BN-609, Inst, for Fluid Dynamics and Appl. Math.,
Univ. of Md., 1969.
[B3] I. Babuska, "Numerical stability in problems of linear algebra,"
SIAM J. Numer. Anal., vol. 9, 1972, pp. 53-77.
[ B4 ] D. W. Bailey and D. E. Crabtree, "Bounds for determinants," Lin.
Alg . ApdI. , vol. 2, 1969, pp . 303-309.
[85] R. E. Bank, "Marching algorithms and block Gaussian elimination,"
in [ B16 , pp. 293-307].
[B6 ] G. H. Barnes, R. M. Brown, M. Kato, D. J. Kuck, D. L. Slotnick,
R. A. Stoker, "The Illiac IV computer," IEEE T-ans . on Como.,
vol. C-17, 1968, pp. 746-757.
[ B 7 ] G. M. Baudet, Asynchronous iterative methods for multiprocessors,
Dept, of Computer Science, Carnegie-Mellon Univ., Nov. 1976.
[B8] F. L. Bauer, "Der Newton- Prozess als quadratisch konvergente
Abkurzung des allgemeinen linearen stationaren Iterationsverfahrens
1. Ordnung (Wittmeyer- Prozess) ," Z. Angew. Math. Mech., vol. 35,
1955, pp. 469-470.
[B9] A. Benson and D. J. Evans, "The successive peripheral block over-
relaxation method," J. Inst. Maths. Applies., vol. 9, 1972,
pp. 68-79.
[BIO] A. Benson and D. J. Evans, "Successive peripheral overrelaxation
and other block methods," J. Comp ■ Phys , , vol. 21, 1976, pp. 1-19.
[Bll] G. Birkhoff and A. George, "Elimination by nested dissection,"
in [Tl, pp. 221-269].
[B12] W. J. Bouknight, S. A. Denenberg, D. E. McIntyre, J. M. Randall,
A. H. Sameh, D. L. Slotnick, "The Illiac IV system," Proc . IEEE ,
vol. 60, 1972, pp. 369-388.
[B13] A. Brandt, Multi-level adaptive techniques I. The multi-grid
method, Rept. RC 6026, IBM T. J. Watson Res. Cen., Yorktown
Heights, N. Y. , 1976.
[B14] J. L. Brenner, "A bound for a determinant with dominant main
diagonal," Proc. AMS, vol. 5, 1954, pp. 631-634.
158
[ B 15 ] C. G. Broyden, "Some condition number bounds for the Gaussian
elimination process," J. Inst. Maths. Applies., vol. 12, 1973,
pp. 273-236.
[B16] J. R. Bunch and D. J. Rose, eds . , Sparse Matrix Computations,
Academic Press, N. Y. , 1976.
[B17] B. L. Buzbee, Application of fast Poisson solvers to the numerical
approximation of parabolic problems. Rent. LA-4950-T, Los Alamos
Sci. Lab., 1972.
[B18] B. L. Buzbee and A. Carasso, "On the numerical computation of para-
bolic problems for preceding times," Math. Comp., vol. 27, 1973,
pp. 237-266.
[B19] B. L. Buzbee, F. W. Dorr, J. A. George and G. H. Golub, "The
direct solution of the discrete Poisson equation on irregular
regions," SIAM J. Numer. Anal., vol. 8, 1971, pp. 722-736.
[B20] B. L. Buzbee, G. H. Golub and C. W. Nielson, "On direct methods
for solving Poisson's equations," S IAM J . Numer . Ana 1 . , vol. 7,
1970, pp. 627-656.
[Cl] A. Carasso, "The abstract backward beam equation," S LAM J . Ma th .
Ana 1 ■ , vol. 2, 1971, pp. '93-212.
[C2] Control Data Corporation, CDC STAR- 100 Instruction execution times,
Arden Hills, Minn., Jan. 1976; in conjunction with measurements
made at Lawrence Livermore Laboratory.
[C3] A. K. Cline, personal communication, 1974.
[C4] P. Concus and G. H. Golub, "Use of fast direct methods for the effi-
cient solution of nonseparable elliptic equations," SLAM J. Numer.
Anal. , vol. 10, 1973, pp. 1103-1120.
[C5] P. Concus, G. H. Golub and D. P. O'Leary, "A generalized conjugate
gradient method for the numerical solution of elliptic partial dif-
ferential equations," in [B16, pp. 309-332].
[C6] S. D. Conte and C. de Boor, Elementary Numerical Analysis, McGraw-
Hill, N. Y., 1972.
M. G. Cox, "Curve fitting with piecewise polynomials," J. Inst.
Maths. Applies., vol. S, 1971, pp. 36-52.
Cray Research, Inc., "Cray- 1 Computer," Chippewa Falls, Wis., 1975.
M. A. Diamond and D. L. V. Ferreira, "On a cyclic reduction method
for the solution of Poisson's equations," SIAM J. Numer. Anal.,
vol. 13, 1976, pp. 54-70.
F. W. Dorr, "The direct solution of the discrete Poisson equation
on a rectangle," SIAM Review, vol. 12, 1970, pp. 248-263.
159
[D3] F. W. Dorr, "An example of ill-conditioning in the numerical solu-
tion of singular perturbation problems," Math. Comp., vol. 25.
1971, pp. 271-283.
[El] S. C. Eisenstat, M. H. Schultz and A. H. Sherman, "Applications of
an element model for Gaussian elimination," in [B16, pp. 85-96].
[E2] D. J. Evans, "A new iterative method for the solution of sparse
systems of linear difference equations," in [R6, pp. 89-100],
[E3] R. K. Even and Y. Wallach, "On the direct solution of Dirichlet's
problem in two dimensions," Computing , vol. 5, 1970, pp. 45-56.
[FI] D. C. Feingold and R. S. Varga, "Block diagonally dominant matrices
and generalizations of the Gerschgorin circle theorem," Pacific J.
Math. , vol. 12, 1962, pp. 1241- 1250.
[F2] D. Fischer, G. Golub, 0. Hald, J. Leiva and 0. Widlund, "On
Fourier-Toeplitz methods for separable elliptic problems," Math ■
Comp . , vol. 28, 1974, pp. 349-368.
[F3] G. von Fuchs, J. R. Roy and E. Schrem, "Hypermatrix solution of
large sets of symmetric positive definite linear equations," Comp .
Math, in Aopl. Mech. and Engng., vol. 1, 1972, pp. 197-216.
[Gl] J. A. George, "Nested dissection of a regular finite element mesh,"
SIAM J. Numer. Anal., vol. 10, 1973, pp. 345-363.
[HI] L. A. Hageman and R. S. Varga, "Block iterative methods for cyclic-
ally reduced matrix equations," Numer. Math., vol. 6, 1964, pp. 106-
119.
[H2] D. E. Heller, "Some aspects of the cyclic reduction algorithm for
block tridiagonal linear systems," SIAM J. Numer. Anal,, vol. 13,
1976, pp. 484-496.
[H3] D. E. Heller, A survey of parallel algorithms in numerical linear
algebra. Dept, of Comp. Sci., Carnegie- Me 11 on Univ. , 1976; SIAM
Review, to appear.
[H4] D. E. Heller, D. K. Stevenson and J. F. Traub, "Accelerated itera-
tive methods for the solution of tridiagonal systems on parallel
computers," J.ACM, vol. 23, 1976, pp. 636-654.
[H5] R. W. Hockney, "A fast direct solution of Poisson's equation using
Fourier analysis," J.ACM, vol. 12, 1965, pp. 95-113.
[H6] R. W. Hockney, "The potential calculation and some applications,"
Methods in Computational Physics, vol. 9, 1970, pp. 135-211.
[H7] E. Horowitz, "The application of symbolic mathematics to a singular
perturbation problem," Proc. ACM Annual Conf., 1972, vol. II,
pp. 816-825.
150
[H8] A. S. Householder, The Theory of Matrices in Numerical Analysis,
Blaisdell, N. Y., 1964.
[Jl] T. L. Jordan, A new parallel algorithm for diagonally dominant
tridiagonal matrices, Los Alamos Sci. Lab., 1974.
[J2] M. L. Juncosa and T. W. Mullikin, "On the increase of convergence
rates of relaxation procedures for elliptic partial difference
equations," J.ACM, vol. 7, 1960, pp. 29-36.
[Kl] W. Kahan, Conserving confluence curbs ill-condition, Dept, of
Comp. Sci., Univ. of Calif., Berkeley, 1972.
[K2] H. T. Kung, "On computing reciprocals of power series, " Numer .
Math. , vol. 22, 1974, pp. 341-348.
[K3] H. T. Kung and J. F. Traub, All algebraic functions can be computed
fast, Dept, of Comp. Sci., Carnegie-Me lion Univ., July 1976.
[LI] J. J. Lambiotte, Jr., The solution of linear systems of equations
on a vector computer, Dissertation, Univ. of Virginia, 1975.
[L2] J. J. Lambiotte, Jr. and R. G. Voigt, "The solution of tridiagonal
linear systems on the CDC STAR- 100 computer," ACM l'rans. Math.
Software, vol. 1, 1975, pp. 308-329.
[Ml] N. Madsen and G. Rodrigue, A comparison of direct methods for tri-
diagonal systems on the CDC STAR-100, Lawrence Livermore Lab.,
1975 .
[M2] M. A. Malcolm and J. Palmer, "A fast direct method for solving a
class of tridiagonal linear systems," CACM, vol. 17, 1974, op. 14-
17.
[M3] W. Miller, "Software for roundoff analysis," ACM Trans. Math.
Software , vol. 1, 1975, pp. 108-128.
[Nl] A. K. Noor and S. J. Voigt, "Hypermatrix scheme for finite element
systems on CDC STAR- 100 computer," Computers and Structures, vol.
5, 1975, pp. 287-296.
[01] J. M. Ortega and R. G. Voigt, Solution of partial differential
equations on vector computers, ICASE, Hamton, Va., 1977.
[02] A. M. Ostrowski, "ijber die Determinanten rait uberwiegender Haupt-
diagonale," Comment. Math. Helv., vol. 10, 1937, pp. 69-96.
[PI] G. Peters and J. H. Wilkinson, "On the stability of Gauss- Jordan
elimination with pivoting," CACM, vol. 18, 1975, pp. 20-24.
[Rl] J. K. Reid, ed., Large Sparse Sets of Linear Equations, Academic
Press, London, 1970.
161
[R2] J. K. Reid, "On the method of conjugate gradients for the solu-
tion of large sparse systems of linear equations," in [Rl,
pp. 231-254].
[R3] F. Robert, "Blocks-H-matrices et convergence des methods itera-
tives classiques par blocks," Lin. Alg . Appl. , vol. 2, 1960,
pp. 223-265. ~
[R4] D. J. Rose and J. R. Bunch, "The role of partitioning in the numer-
ical solution of sparse systems," in [R6, pp. 177-187],
[R5] D. J. Rose and G. F. Whitten, "A recursive analysis of dissection
strategies," in [B16, pp. 59-83].
[R6] D. J. Rose and R. A. Willoughby, eds . , Sparse Matrices and Their
Applications , Plenum Press, N. Y., 1972.
[R7] J. B. Rosser, The direct solution of difference analogs of
Poisson's equation, rep. MRC-TSR- 797 , Math. Res. Cen., Madison,
Wis., 1967.
[51] J. Schroder and U. Trottenberg, "Reduktionsverfahren fur Differ-
enzengleichungen bei Randwertaufgaben I," N'uner . Math., vol. 22,
1973, pp. 37-68.
[Sla] J. Schroder, U. Trottenberg and H. Reutersberg, "Reduktionsver-
fahren fur Dif f erenzengleichungen bei Randwertaufgaben II," Numer .
Ma th . , vol. 16, 1976, pp. 429-459.
[52] H. S. Stone, "ParalleL tridiagonal equation solvers," ACM Trans ■
Math. Software, vol. 1, 1975, pp. 289-307.
[53] H. S. Stone, ed., Introduction to Computer Architecture, Science
Research Associates, Palo Alto, Calif., 1975.
[54] G. Strang and G. J. Fix, An Analysis of the Finite Element Method,
Prentice-Hall, Englewood Cliffs, N. J. , 1973.
[55] P. N. Swarztrauber , "The direct solution of the discrete Poisson
equation on the surface of a sphere," J. Comp. Phys., vol. 15,
1974, pp. 46-54.
[56] P. N. Swarztrauber, "A direct method for the discrete solution of
separable elliptic equations," SIAM J. Numer. Anal., vol. 11,
1974, pp. 1136-1150.
[57] P. N. Swarztrauber and R. A. Sweet, "The direct solution of the
discrete Poisson equation on a disc," SIAM J. Numer. Anal., vol.
10, 1973, pp. 900-907.
[58] P. Swarztrauber and R. Sweet, Efficient FORTRAN subprograms for
the solution of elliptic partial differential equations, Rept.
NCAR-TN/lA- 109, Nat. Cen. for Atmospheric Res., Boulder, Col.,
1975.
162
[S9] R. A. Sweet, "Direct methods for the solution of Poisson's equa-
■ tion on a staggered grid," J . Comp ■ Phys . , vol. 12, 1973, do. 422-
42S .
[510] R. A. Sweet, "A generalized cyclic reduction algorithm," SIAM J.
Numer ■ Anal . , vol. 11, 1974, pp. 506-320.
[ 5 1 1 ] R. A. Sweet, "A cyclic reduction algorithm for block tridiagonal
systems of arbitrary dimension," contributed paper, Symp . on
Sparse Matrix Computations, Argonne Nat. Lab., Sept. 1975.
[Tl] J. F. Traub, ed.. Complexity of Sequential and Parallel Numerical
Algorithms , Academic Press, N. Y. , 1973.
[T2] J. F. Traub, "Iterative solution of tridiagonal systems on parallel
or vector computers," in [Tl, pp. 49-82].
[T3] J. F. Traub, ed. , Analytic Computational Complexity, Academic
Press, N. Y., 1976.
[VI] J. M. Varah, "On the solution of block tridiagonal systems arising
from certain finite difference equations." Math. Comp., vol. 26,
1972, pp. 359-868.
[V2] J. M. Varah, "A comparison of some numerical methods for two-point
boundary value problems,” Math . Comp ■ , vol. 28, 1974, pp. 743-755.
[V3] R. S. Varga, Matrix Iterative Analysis, Prentice-Hall, Englewood
Cliffs, N. J., 1962.
[V4] R. S. Varga, "On recurring theorems on diagonal dominance," Lin .
Alg . Appl ■ , vol. 13, 1976, pp. 1-9.
[Wl] W. J. Watson, "The Tl ASC, a highly modular and flexible super-
computer architecture," AFIPS Fall 1972, AFIPS Press, Montvale,
N. J., vol. 41, pt. 1, pp. 221-229.
[W2 ] 0. B. Widlund, "On the use of fast methods for separable finite
difference equations for the solution of general elliptic problems,"
in [R6 , pp. 121-131].
[W3 ] J. H. Wilkinson, "Error analysis of direct methods of matrix inver-
sion," J.ACM, vol. 8, 1961, pp. 281-330.
[W4 ] W. A. Wulf and C. G. Bell, "C.mmp, a multi-mini-processor," AFIPS
Fall 1972, AFIPS Press, Montvale, N. J., vol. 41, pt. 2, pp. 765-
777 .
[yl] D. Y. Y. Yun, "Hensel meets Newton - algebraic constructions in an
analytic setting," in [T3, pp. 205-215].
II
SECURITY CLASSIFICATION OF THIS PAGE (Whtn Data Enter'd)
REPORT DOCUMENTATION PAGE
t. REPORT NUMBER
4 TITLE (and Su6f/f/e)
READ INSTRUCTIONS
BEFORE COMPLETING FORM
3. RECIPIENT'S CATALOG NUMBER
5. TYPE OF REPORT & PERIOD COVERED
DIRECT AND ITERATIVE METHODS FOR BLOCK
TRIDIAGONAL LINEAR SYSTEMS
Interim
6 PERFORMING ORG. REPORT NUMBER
7. AUTHORS.)
Don Eric Heller
9- PERFORMING ORGANIZATION NAME AND ADDRESS
Carnegie-Mellon University
Computer Science Dept.
i s?i
II. CONTROLLING OFFICE NAME AND ADDRESS
Office of Naval Research
Arlington, VA 22217
8. CONTRACT OR GRANT NUMBERf*.)
MCS75-222-55
N00014- 76-C-0370 ;
NR 044-422
10 PROGRAM ELEMENT. PROJECT, TASK
AREA 4 WORK UNIT NUMBERS
12 REPORT DATE
April 1977
13, NUMBER OF PAGES
169
14 MONITORING AGENCY NAME 9 ADDRESS (II dltlerent front Controlling Office) IS SECURITY CLASS, fo I t hie report)
UNCLASSIFIED
15a, DECLASSIFICATION DOWNGRADING
SCHEDULE
16. DISTRIBUTION STATEMENT (of this Report)
Approved for public release; distribution unlimited.
17 DISTRIBUTION STATEMENT (of the abstract entered In Block 20, It different from Report)
19 KEY WORDS (Continue on reverse side If necessary and Identify by block number)
20 ABSTRACT (Continue on reverse side It necessary and Identity by block number) Block tridiagOnal Systems of
linear equations occur frequently in scientific computations, often forming the
:ore of more complicated problems. Numerical methods for solution of such sys-
:ems are studied with emphasis on efficient methods for a vector computer. A
:onvergence theory for direct methods under conditions of block diagonal domi-
lace is developed, demonstrating stability, convergence and approximation prop-
irties of direct methods. Block elimination (LU factorization) is linear,
:yclic odd-even reduction is quadratic, and higher-order methods exist. The
dd-even methods are variations of the Quadratic Newton iteration for the
DO , J AN^73 1473 edition of i mov «5 is obsolete continued
S/N 0102-014- 6601 ‘
SECURITY CLASSIFICATION of THIS PAGE f«7i»n Dele Entered)
UNCLASSIFIED
liuHITV CLASSIFICATION or THIS PAGEfl »hen Palm Enfrid)
20. abstract continued
inverse matrix, and are the only quadratic methods within a certain reasonable
class of algorithms. Semi-direct methods based on the quadratic convergence
of odd-even reduction prove useful in combination with linear iterations for
an approximate solution. An execution time analysis for a pipeline computer
is given, with attention to storage requirements and the effect of machine
constraints on vector operations.
UNCLASSIFIED
SECURITY CLASSIFICATION OF THIS PAGEnr?i»n Oafs Enfrtd)
Live Music Archive
Librivox Free Audio
Metropolitan Museum
Cleveland Museum of Art
Internet Arcade
Console Living Room
Books to Borrow
Open Library
TV News
Understanding 9/11